Natural cubic spline interpolation - check and suggest better way

945 Views Asked by At

I was given the following interpolation nodes:

$(0,10),(\frac{1}{2},8),(1,5),(2,2),(3,1)$

and I was asked to find the natural cubic spline interpolation between every 2 points.

I want to show you what I did, and hopefully someone will suggest a better way, since this is very long and very tedious, and I only managed to accomplish $3$ points. I will show you.

Suppose in $[0,\frac{1}{2}]$ our interpolating spline is $\xi(x)=ax^3+bx^2+cx+d$

notice that $\xi(x)'=3ax^2+2bx+c$ and that $\xi(x)''=6ax+2b$.

Since we were asked for a natural cubic spline, and $x_0=0$, we will demand $\xi(0)''=2b=0$ and so $b=0$.

we know that $\xi(0)=d=10$ and so $d=10$, and we also know that $\xi(\frac{1}{2})=\frac{1}{8}a+\frac{1}{4}b+\frac{1}{2}c+d=\frac{1}{8}a+\frac{1}{2}c+10=8$

so we can infer that $a=-16-4c$

So all in all after we use all the information, our first cubic polynomial is $\xi(x)=(-16-4c)x^3+cx+10$ where $x \in [0,\frac{1}{2}]$

its derivative is $\xi(x)'=(-48-12c)x^2+c$

And the second derivative is $\xi(x)''=(-96-24c)x$

Now let's look at $(\frac{1}{2},1]$: Suppose our cubic polynomial there is

$\psi(x)=ex^3+fx^2+gx+h$

Notice that $\psi(x)'=3ex^2+2fx+g$ and that $\psi(x)''=6ex+2f$

since we demand that the second derivative is continuous at the interpolation nodes, we must demand

$\psi(\frac{1}{2})''=(-96-24c)\frac{1}{2}=-48-12c=\psi(\frac{1}{2})''$

So $\psi(\frac{1}{2})''=3e+2f=-48-12c$ and so $f=-24-6c-\frac{3}{2}e$.

This insures that the second derivative is continuous, but in addition to that, we must demand that the first derivative is continuous at the interpolation nodes. and so

$\xi(\frac{1}{2})'=(-48-12c)\frac{1}{4}+c=-12-3c+c=-12-2c=\psi(\frac{1}{2})'$, and so

$\psi(\frac{1}{2})'=\frac{3}{4}e+f+g=\frac{3}{4}e-24-6c-\frac{3}{2}e+g=-\frac{3}{4}e-24-6c+g=-12-2c$

and so we can infer $g=\frac{3}{4}e+12+4c$

So right now what we now of $\psi$ is that

$\psi(x)=ex^3+(-24-6c-\frac{3}{2}e)x^2+(\frac{3}{4}e+12+4c)x+h$

However, we also know that $\psi(\frac{1}{2})=8$, and so

$\psi(\frac{1}{2})=\frac{e}{8}-6-\frac{3}{2}c-\frac{3}{8}e+\frac{3}{8}e+6+2c+h=\frac{e}{8}+\frac{1}{2}c+h=8$

And so we can infer

$h=8-\frac{e}{8}-\frac{1}{2}c$

So right now $\psi(x)=ex^3+(-24-6c-\frac{3}{2}e)x^2+(\frac{3}{4}e+12+4c)x+8-\frac{e}{8}-\frac{1}{2}c$

However, we also know that $\psi(1)=5$, and so

$\psi(1)=e-24-6c-\frac{3}{2}e+\frac{3}{4}e+12+4c+8-\frac{e}{8}-\frac{1}{2}c=\frac{1}{8}e-4-\frac{5}{2}c=5$

And from here we can infer $e=72+20c$

Now if we only want to work with $c$:

$f=-24-6c-\frac{3}{2}e=-24-6c-108-30c=-132-36c$

$g=\frac{3}{4}e+12+4c=54+15c+12+4c=66+19c$

$h=8-\frac{e}{8}-\frac{1}{2}c=8-9-\frac{5}{2}c-\frac{1}{2}c=-1-3c$

And so overall

$\xi(x) =(-16-4c)x^3+cx+10$ where $x \in [0,\frac{1}{2}]$

And

$\psi(x)=(72+20c)x^3+(-132-36c)x^2+(66+19c)x-1-3c$ where $x \in (\frac{1}{2},1]$

Now, after all that is done, I have a few questions.

1) Can I choose any $c$ that I like? Will there always be a "free" coefficient?

2) Is there a better and faster way to do this? This was very tiring, and we only took care of half the points! There is still $(1,2]$ and $(2,3]$ left!

Note: The answer in my notebook is different, and is correct. This means that somewhere on this page I was wrong. this only strengthens my point that this method is tiring, long, and very tedious, which is fertile soil for silly mistakes.

1

There are 1 best solutions below

0
On BEST ANSWER

Here is an algorithm for points $(x_0,y_0),(x_1,y_1),\ldots,(x_n,y_n)$:

  1. Calculate for $i\in[0,n-1]$:

    • $\displaystyle w_{i} = x_{i+1}-x_{i}$

    • $\displaystyle h_{i} = \frac{y_{i+1}-y_{i}}{w_{i}}$

  2. Define:

    • $\displaystyle f''_{0} = 0$

    • $\displaystyle f''_{n} = 0$

  3. Calculate for $i\in[0,n-2]$:

    • $\displaystyle f''_{i+1} = 3\frac{h_{i+1}-h_{i}}{w_{i+1}+w_{i}}$
  4. Calculate for $i\in[0,n-1]$:

    • $\displaystyle A_{i} = \frac{f''_{i+1}-f''_{i}}{6w_{i}}$

    • $\displaystyle B_{i} = \frac{f''_{i}}{2}$

    • $\displaystyle C_{i} = h_{i}-w_{i}\frac{f''_{i+1}+2f''_{i}}{6}$

    • $\displaystyle D_{i} = y_{i}$

Please note that the algorithm above assumes $x_0 < x_1 < \ldots < x_n$.