Numerical Computation Interpolation

Data and interpolating functions

DEFINITION: The function $y=P(x)$ interpolates the data points $(x_1,y_1),\dots ,(x_n,y_n)$ if $P(x_i)=y_i$, for each $1\le i\le n.$

In polynomial interpolation, we find a polynomial (or set of polynomials) that passes through each data point.

Lagrange Interpolation

Given two (different) data points $(x_0,y_0)$ and $(x_1,y_1)$, there is a unique line through these points. The equation of the line through $(x_0,y_0)$ and $(x_1,y_1)$ is

$$y-y_0=\frac{y_1-y_0}{x_1-x_0}(x-x_0).$$

NOTE: Consider the function

$$p(x)=y_1\frac{x-x_0}{x_1-x_0}+y_0\frac{x-x_1}{x_0-x_1}.$$

Clearly, $p(x_0)=y_0$ and $p(x_1)=y_1$. Also, $p(x)$ is linear (degree = 1). Thus, $p(x)$ must also be the equation of the line through $(x_0,y_0)$ and $(x_1,y_1)$. The terms in $p(x)$ are called Lagrange Polynomials.

DEFINITION: For a given set of $n+1$ nodes $x_i$, the Lagrange polynomials are the $n+1$ polynomials $L_i$ defined by

$$L_i(x)=\prod _{j=0,j\ne i}^n\frac{x-x_j}{x_i-x_j}$$

DEFINITION: We define a Lagrange Interpolating polynomial $p_n(x)$ by

$$p_n(x)=\sum _{i=0}^n{y}_i{L}_i(x)$$

EXAMPLE

Find an interpolating polynomial for the data points $(0,1)$, $(2,2)$ and $(3,4)$.

SOLUTION:

$$\begin{array}{rl}{P}_2(x)& =1\frac{(x-2)(x-3)}{(0-2)(0-3)}+2\frac{(x-0)(x-3)}{(2-0)(2-3)}+4\frac{(x-0)(x-2)}{(3-0)(3-2)}\\ & =\frac{1}{6}(x^2-5x+6)+2(-\frac{1}{2})(x^2-3x)+4\left(\frac{1}{3}\right)(x^2-2x)\\ & =\frac{1}{2}{x}^2-\frac{1}{2}x+1.\end{array}$$

THEOREM:
Let $(x_1,y_1),\dots ,(x_n,y_n)$ be $n$ points in the plane with distinct $x_i$. Then there exists one and only one polynomial $P$ of degree $n-1$ or less that satisfies $P(x_i)=y_i$ for $i=1,\dots ,n$.

PROOF:
The existence is proved by the explicit formula for Lagrange interpolation.

To show there is only one, assume for the sake of argument that there are two, say, $P(x)$ and $Q(x)$, that have degree at most $n-1$ and that both interpolate all $n$ points. That is, we are assuming that $P(x_1)=Q(x_1)=y_1$, $P(x_2)=Q(x_2)=y_2$, $\dots$, $P(x_n)=Q(x_n)=y_n$. Now define the new polynomial $H(x)=P(x)-Q(x)$. Clearly, the degree of $H$ is also at most $n-1$, and note that $0=H(x_1)=H(x_2)=\cdots =H(x_n)$; that is, $H$ has $n$ distinct zeros (roots). According to the Fundamental Theorem of Algebra, a degree $d$ polynomial can have at most $d$ zeros, unless it is the identically zero polynomial. Therefore, $H$ is the identically zero polynomial, and $P(x)\equiv Q(x)$. We conclude that there is a unique $P(x)$ of degree $\le n-1$ interpolating the $n$ points $(x_i,y_i)$.

Newton Interpolation

In Newton’s method, we create a polynomial $p(x)$ for data $(x_i,y_i)$ such that $p(x_i)=y_i$ (Just as we did in the Lagrange Method).
The difference is that we construct $p(x)$ iteratively.
That is we construct a sequence of polynomials $p_0(x),p_1(x),\dots ,p_n(x)$ such that $p_n(x)$ is the desired interpolating polynomial.

Algorithm:

1. Let $p_0(x)=y_0$ ($p_0(x)$ matches data at $x_0$)
2. Suppose we have calculated $p_k(x)$ then, set$$p_{k+1}(x)=p_k(x)+c(x-x_0)(x-x_1)\dots (x-x_k)$$for some constant $c$.

NOTE:
$p_{k+1}(x_i)=p_k(x_i)$ for $i=0,1,\dots ,k$. So, $p_{k+1}(x_i)$ correctly interpolates all data up to $x_k$. To make it match at $x_{k+1}$ we set

$$y_{k+1}=p_k(x_{k+1})+c(x_{k+1}-x_0)(x_{k+1}-x_1)\dots (x_{k+1}-x_k)$$

and solve for $c$.

EXAMPLE
Find an interpolating polynomial for the data points $(0,1)$, $(2,2)$ and $(3,4)$ using Newton interpolating method.

$$p_1(x)=c_1$$

Let $1=p_1(0)=c_1\Rightarrow c_1=1$

$$p_2(x)=p_1(x)+c_2(x-0)=1+c_2(x-0)$$

Let $2=p_2(2)=1+2c_2\Rightarrow c_2=\frac{1}{2}$

$$p_3(x)=p_2(x)+c_3(x-0)(x-2)=1+\frac{1}{2}x+c_{3}x(x-2)$$

Let $4=p_3(3)=1+\frac{3}{2}+3c_3\Rightarrow c_3=\frac{1}{2}$

Hence,

$$p(x)=p_3(x)=\frac{1}{2}{x}^2-\frac{1}{2}x+1$$

which is the same as the result we found by Lagrange interpolating method.

Newton’s Nested formula

Newton’s Method vs Lagrange:

• Both methods produce the same polynomial.
• Newton’s Method is more flexible. It is easy to add new data points and get a new interpolating polynomial.
• The coefficients in Newton’s Method (determined by solving for the constants $c$) can be calculated very efficiently. (Newton divide difference)

From previous example, without calculating the coefficient, we can write

$$p(x)=c_1+c_2(x-0)+c_3(x-0)(x-2).$$

It can be rewritten as

$$p(x)=c_1+(x-0)[c_2+c_3(x-2)].$$

In general,

$$p_n(x)=c_1+[c_2(x-x_1)]+\cdots +[c_n(x-x_1)(x-x_1)\dots (x-x_{n-1})]$$

can be rewritten as

$$p_n(x)=c_1+(x-x_1)[c_1+(x-x_2)[c_3+(x-x_3)[\dots ]]]$$

called Newton’s nested formula.

$$p_n(x)=c_1+(x-x_1)[c_2+(x-x_2)[c_3+(x-x_3)[\dots ]]]$$

allows a better way of calculating $p_n(t)$ for a value $t$ (or $x$):

$$\begin{array}{rl}{v}_1& =c_n\\ v_2& =c_{n-1}+(t-x_{n-1})v_1\\ v_3& =c_{n-2}+(t-x_{n-2})v_2\\ & ⋮\\ v_n& =c_1+(t-x_1)v_{n-1}\end{array}$$

Then, $v_n=p_n(t)$ Newton’s Nested formula requires $2(n-1)$ additions and $n-1$ multiplications. $3(n-1)\sim O(n)$ Lagrange Interpolation requires at least $n^2$ operations $\sim O(n^2)$

$$p(x)=\sum _k\left(\prod _{j\ne k}^n\frac{x-x_j}{x_k-x_j}\right)y_k$$

Following we will introduce an efficient way to calculate the coefficient $c_k$.

DEFINITION:

For a given collection of data $\{(x_i,f(x_i))\}$ a $k$th order divided difference is a function of $k+1$ (not necessarily distinct) data values written as $f[x_i,x_{i+1},...,x_{i+k}]$.

$1^{st}$ Order: $f[x_i]=f(x_i)=y_i$$k^{th}$ Order:

$$f[x_i,x_{i+1},...,x_{i+k-1},x_{i+k}]=\frac{f[x_{i+1},x_{i+2},...,x_{i+k}]-f[x_i,x_{i+1},...,x_{i+k-1}]}{x_{i+k}-x_i}$$

THEOREM

For Newton’s Interpolation formula

$$p_n(x)=c_1+[c_2(x-x_1)]+\cdots +[c_n(x-x_1)(x-x_1)\dots (x-x_{n-1})]$$

the coefficients $c_k$ can be calculated as

$$c_k=f[x_1,x_2,...,x_k]$$

The recursive definition of the Newton’s divided differences allows arrangement into a convenient table. For three points the table has the form

$$\begin{array}{cccc}{x}_1& \stackrel{c_1}{f[x_1]}& & \\ & & \stackrel{c_2}{f[x_1{x}_2]}& \\ x_2& f[x_2]& & \stackrel{c_3}{f[x_1{x}_2{x}_3]}\\ & & f[x_2{x}_3]& \\ x_3& f[x_3]& & \end{array}$$

The coefficients of the Newton’s Interpolation formula can be read from the top edge of the triangle.

EXAMPLE

Use divided differences to find the interpolating polynomial passing through the points $(0,1)$, $(2,2)$ and $(3,4)$.

SOLUTION

The divided differences are calculated as follows:

$$\frac{2-1}{2-0}=\frac{1}{2}$$$$\frac{4-2}{3-2}=2$$$$\frac{2-\frac{1}{2}}{3-0}=\frac{1}{2}$$

The coefficients of the polynomial $1$, $1/2$, $1/2$ can be read from the top edge of the table. The interpolating polynomial can be written as

$$p(x)=1+\frac{1}{2}(x-0)+\frac{1}{2}(x-0)(x-2),$$

or, in nested form

$$p(x)=1+(x-0)(\frac{1}{2}+\frac{1}{2}(x-2)).$$

Approximating Polynomials and Interpolation Error

A major use of polynomial interpolation is to replace evaluation of a complicated function by evaluation of a polynomial, which involves only elementary computer operations like addition, subtraction, and multiplication.

EXAMPLE:

Interpolate the function $f(x)=\sin x$ at 4 equally spaced points on $[0,\pi /2]$.

SOLUTION

$$p(x)=0+x(0.9549+(x-\frac{\pi }{6})(-0.2443+(x-\frac{\pi }{3})(-0.1139))).$$$$\begin{array}{rllll}0& 0.0000& & & \\ & & 0.9549& & \\ \pi /6& 0.5000& & -0.2443& \\ & & 0.6990& & -0.1139\\ 2\pi /6& 0.8660& & -0.4232& \\ & & 0.2559& & \\ 3\pi /6& 1.0000& & & \end{array}$$

We can repeat the process on the interval $[\frac{\pi }{2},\pi ],[\pi ,\frac{3\pi }{2}],...$ based on $\sin x$ symmetry and antisymmetric.

The error for degree 3 interpolation of $\sin x$ is shown below.

In order to get enough correct digits to fill the calculator readout, we’ll need to know a little more about interpolation error.

[A graph is shown here with the function $y=\sin (x)$ and its interpolation, along with a table of values and errors]

$x$	$\sin x$	$\text{sin1}(x)$	error
1	0.8415	0.8411	0.0004
2	0.9093	0.9102	0.0009
3	0.1411	0.1428	0.0017
4	$-0.7568$	$-0.7557$	0.0011
14	0.9906	0.9928	0.0022
1000	0.8269	0.8263	0.0006

THEOREM

Assume that $P(x)$ is the (degree $n-1$ or less) interpolating polynomial fitting the $n$ points $(x_1,y_1),...,(x_n,y_n)$. The interpolation error is

$$f(x)-P(x)=\frac{(x-x_1)(x-x_2)...(x-x_n)}{n!}{f}^{(n)}(c)$$

where $c$ lies between the smallest and largest of the number $x,x_1,...,x_n$.

PROOF

Consider adding one more point $x$ to the set of interpolation points. The new interpolation polynomial would be

$$P_n(t)=P_{n-1}(t)+f[x_1,x_2,...,x_n,x](t-x_1)\cdots (t-x_n).$$

At the point $x$, $P_n(x)=f(x)$, so

$$f(x)=P_{n-1}(x)+f[x_1,x_2,...,x_n,x](x-x_1)\cdots (x-x_n).$$

Now define

$$h(t)=f(t)-P_{n-1}(t)-f[x_1,x_2,...,x_n,x](t-x_1)\cdots (t-x_n).$$

Note that $h(x)=0=h(x_1)=\cdots =h(x_n)$.

By Rolle’s theorem, there are $n$ points between each neighboring pair of $n+1$ points $x_1,x_2,...,x_n,x$ that $h^{\prime }=0$.

Furthermore, there are $n-1$ points between each neighboring pair of those $n$ points found above that $h^{″}=0$.

Continuing in this way, there must be one point $c$ for which $h^{(n)}=0$, where $c$ lies between the smallest and largest of the number $x,x_1,...,x_n$.

Takes $n$th order derivative for above equation, we got

$$h^{(n)}(t)=f^{(n)}(t)-n!f[x_1,x_2,...,x_n,x].$$

because the $n$th order derivative of $P_{n-1}(t)$ is zero. Substituting $c$ gives

$$f[x_1,x_2,...,x_n,x]=\frac{f^{(n)}(c)}{n!},$$

which leads to

$$f(x)=P_{n-1}(x)+\frac{f^{(n)}(c)}{n!}(x-x_1)\cdots (x-x_n).$$

EXAMPLE

Find an upper bound for the difference at $x=0.25$ and $x=0.75$ between $f(x)=e^x$ and the polynomial that interpolates it at the points $-1,-0.5,0,0.5,1$.

SOLUTION

$$f(x)-P_4(x)=\frac{(x+1)(x+\frac{1}{2})x(x-\frac{1}{2})(x-1)}{5!}{f}^{(5)}(c),$$

which $-1<c<1$. So $f^{(5)}(c)=e^c\le e$.

At $x=0.25$,

$$|e^{0.25}-P_4(0.25)|\le \left|\frac{(1.25)(0.75)(0.25)(-0.25)(-0.75)}{5!}e\right|\approx 0.000995.$$

At $x=0.75$,

$$|e^{0.75}-P_4(0.75)|\le \left|\frac{(1.75)(1.25)(0.75)(0.25)(-0.25)}{5!}e\right|\approx 0.002323.$$

Cubic Splines

A spline curve is a curve that is made of a set of simple curves (lines, quadratics, cubic) that are joined together. A spline is piece-wise defined. That is, it is defined over a set of sub-intervals of a given interval. This set is called a partition.

DEFINITION: A partition of an interval $[a,b]$ is a sequence $\{t_i{\}}_{i=1}^n$ of points between $a$ and $b$ such that

$$a=t_1<t_2<\cdots <t_n=b$$

The numbers $t_i$ ($i=2$ to $n-1$) are called knots.

DEFINITION: A function $S$ is a linear spline on $[a,b]$ if

• The domain of $S$ is $[a,b]$
• $S$ is continuous on $[a,b]$
• There is a partition $\{t_i{\}}_{i=1}^n$ of points on $[a,b]$ such that $S$ is a linear function on each sub-interval $[t_i,t_{i+1}]$

---

A linear spline is defined by its values at the set of knots. Given the table of values there is a unique linear spline with those values.

$t$	$t_1$	$t_2$	$\cdots$	$t_n$
$y$	$y_1$	$y_2$	$\cdots$	$y_n$

On each sub-interval, $[t_i,t_{i+1}]$, the linear spline is defined by

$$S_i(x)=y_i+\frac{y_{i+1}-y_i}{t_{i+1}-t_i}(x-t_i)$$

Then, $S_i(t_i)=y_i$ and $S_i(t_{i+1})=y_{i+1}$

EXAMPLE:

$t$	0.0	0.1	0.4	0.5	0.75	1.0
$y$	1.3	4.5	2.0	2.1	5.0	3.0

[A graph of a linear spline connecting the points from the table is shown here]

---

Quadratic Spline

DEFINITION: A function $Q$ is a quadratic spline on $[a,b]$ if

• The domain of $Q$ is $[a,b]$
• $Q$ is continuous on $[a,b]$
• $Q^{\prime }$ is continuous on $[a,b]$
• There is a partition $\{t_i{\}}_{i=1}^n$ of points on $[a,b]$ such that $Q$ is a polynomial of degree $\le 2$ on each sub-interval $[t_i,t_{i+1}]$

$$\begin{array}{rlr}Q(x)& =a_1+b_1(x-t_1)+c_1(x-t_1{)}^2,& t_1\le x\le t_2\\ & =a_2+b_2(x-t_2)+c_2(x-t_2{)}^2,& t_2\le x\le t_3\\ & =⋮\\ & =a_{n-1}+b_{n-1}(x-t_{n-1})+c_{n-1}(x-t_{n-1}{)}^2,& t_{n-1}\le x\le t_n\end{array}$$

Thus, we need to determine $3(n-1)$ different constants.

---

We know that:

$$\begin{array}{rl}{Q}_i(t_i)& =y_i\\ Q_i(t_{i+1})& =y_{i+1}\\ Q_i^{\prime }(t_{i+1})& =Q_{i+1}^{\prime }(t_{i+1})(Q^{\prime }\text{ is continuous at knots})\end{array}$$

This gives $(n-1)+(n-1)+(n-2)=3(n-1)-1$ equations for $3(n-1)$ unknown!

Need one more condition on $Q$. Could use one of these

$$\begin{array}{rl}{c}_1& =0(Q_1\text{ is a line})\\ Q^{\prime }(t_1)& =z(\text{slope at start})\\ Q^{″}(t_1)& =w(\text{curvature at start})\end{array}$$

---

EXAMPLE Data $x=[014]$, $y=[1-21]$. Find $Q_1(x)$ and $Q_2(x)$ when $Q_1^{\prime }(0)=0$. SOLUTION Define

$$\begin{array}{rlr}{Q}_1(x)& =a_1+b_1(x-0)+c_1(x-0{)}^2,& 0\le x\le 1\\ Q_2(x)& =a_2+b_2(x-1)+c_2(x-1{)}^2,& 1\le x\le 4\\ Q_1^{\prime }(x)& =b_1+2c_{1}x\\ Q_2^{\prime }(x)& =b_2+2c_2(x-1)\end{array}$$

By quadratic spline definition,

$$\begin{array}{rl}& Q_1^{\prime }(0)=0\Rightarrow b_1=0;Q_1(0)=1\Rightarrow a_1=1;Q_2(1)=-2\Rightarrow a_2=-2\\ & Q_1(1)=-2\Rightarrow c_1+1=-2\Rightarrow c_1=-3\\ & Q_1^{\prime }(1)=Q_2^{\prime }(1)\Rightarrow b_2=-6\\ & Q_2(4)=1\Rightarrow 9c_2+3b_2+a_2=1\Rightarrow c_2=\frac{7}{3}.\end{array}$$

Hence

$$\begin{array}{rl}{Q}_1(x)& =-3x^2+1\\ Q_2(x)& =\frac{7}{3}(x-1{)}^2-6(x-1)-2\end{array}$$

---

• 计算步长： ${\delta }_i=x_{i+1}-x_i$。
• 利用插值条件： 直接得出 $a_i=y_i$。
• 利用连续性： 列出关于 $b_i$ 和 $c_i$ 的递归关系。
• 利用一阶导数连续： 得到 $b_{i+1}=b_i+2c_i{\delta }_i$。
• 代入边界条件： 假设 $c_0=0$，然后从左到右依次解出所有的 $b_i$ 和 $c_i$。

---

Cubic Spline

DEFINITION: A function $S$ is a cubic spline on $[a,b]$ if

• The domain of $S$ is $[a,b]$
• $S$ is continuous on $[a,b]$
• $S^{\prime }$ is continuous on $[a,b]$
• $S^{″}$ is continuous on $[a,b]$
• There is a partition $\{t_i{\}}_{i=1}^n$ of points on $[a,b]$ such that $S$ is a polynomial of degree $\le 3$ on each sub-interval $[t_i,t_{i+1}]$

$$\begin{array}{rlr}S(x)& =y_1+b_1(x-t_1)+c_1(x-t_1{)}^2+d_1(x-t_1{)}^3,& t_1\le x\le t_2\\ & =y_2+b_2(x-t_2)+c_2(x-t_2{)}^2+d_2(x-t_2{)}^3,& t_2\le x\le t_3\\ & =⋮\\ & =y_{n-1}+b_{n-1}(x-t_{n-1})+c_{n-1}(x-t_{n-1}{)}^2+d_{n-1}(x-t_{n-1}{)}^3,& t_{n-1}\le x\le t_n\end{array}$$

---

Property 1:$S_i(t_i)=y_i$ and $S_i(t_{i+1})=y_{i+1}$ for $i=1,...,n-1$.

Property 2:$S_{i-1}^{\prime }(t_i)=S_i^{\prime }(t_i)$ for $i=2,...,n-1$.

Property 3:$S_{i-1}^{″}(t_i)=S_i^{″}(t_i)$ for $i=2,...,n-1$.

---

EXERCISE Check that

$$S_1(x)=2-\frac{13}{8}(x-1)+0(x-1{)}^2+\frac{5}{8}(x-1{)}^3\text{ on }[1,2]$$$$S_2(x)=1+\frac{1}{4}(x-2)+\frac{15}{8}(x-2{)}^2-\frac{5}{8}(x-2{)}^3\text{ on }[2,4]$$$$S_3(x)=4+\frac{1}{4}(x-4)-\frac{15}{8}(x-4{)}^2+\frac{5}{8}(x-4{)}^3\text{ on }[4,5]$$

satisfies all cubic spline properties for the data points $(1,2),(2,1),(4,4)$ and $(5,3)$.

---

Constructing a spline from a set of data points means finding the coefficients $b_i,c_i,d_i$ that make Properties 1-3 hold. Properties 1-3 gives $n-1+2(n-2)=3n-5$ independent equations. The total number of the coefficients $b_i,c_i,d_i$ is $3(n-1)=3n-3$. Therefore, solving for the coefficients is a problem of solving $3n-5$ linear equations in $3n-3$ unknowns. Two extra equations needed!

Property 4 Natural spline$S_1^{″}(t_1)=0$ and $S_{n-1}^{″}(t_n)=0$.

---

Property 1: $S_i(t_{i+1})=y_{i+1}$ for $i=1,...,n-1$. (n-1 equations)

$$\begin{array}{rl}{y}_2& =S_1(t_2)=y_1+b_1(t_2-t_1)+c_1(t_2-t_1{)}^2+d_1(t_2-t_1{)}^3\\ & ⋮\\ y_n& =S_{n-1}(t_n)=y_{n-1}+b_{n-1}(t_n-t_{n-1})+c_{n-1}(t_n-t_{n-1}{)}^2+d_{n-1}(t_n-t_{n-1}{)}^3.\end{array}$$

Property 2: $S_{i-1}^{\prime }(t_i)=S_i^{\prime }(t_i)$ for $i=2,...,n-1$. (n-2 equations)

$$\begin{array}{rl}0=S_1^{\prime }(t_2)-S_2^{\prime }(t_2)& =b_1+2c_1(t_2-t_1)+3d_1(t_2-t_1{)}^2-b_2\\ & ⋮\\ 0=S_{n-2}^{\prime }(t_{n-1})-S_{n-1}^{\prime }(t_{n-1})& =b_{n-2}+2c_{n-2}(t_{n-1}-t_{n-2})+3d_{n-2}(t_{n-1}-t_{n-2}{)}^2-b_{n-1}.\end{array}$$

Property 3: $S_{i-1}^{″}(t_i)=S_i^{″}(t_i)$ for $i=2,...,n-1$. (n-2 equations)

$$\begin{array}{rl}0=S_1^{″}(t_2)-S_2^{″}(t_2)& =2c_1+6d_1(t_2-t_1)-2c_2\\ & ⋮\\ 0=S_{n-2}^{″}(t_{n-1})-S_{n-1}^{″}(t_{n-1})& =2c_{n-2}+6d_{n-2}(t_{n-1}-t_{n-2})-2c_{n-1}.\end{array}$$

---

Let ${\delta }_i=t_{i+1}-t_i$ and ${\mathrm{\Delta }}_i=y_{i+1}-y_i$. If we know all $c_i$, from property 3 we can solve

$$d_i=\frac{c_{i+1}-c_i}{3{\delta }_i}$$

From property 1 we can solve

$$\begin{array}{rl}{b}_i& =\frac{{\mathrm{\Delta }}_i}{{\delta }_i}-c_i{\delta }_i-d_i{\delta }_i^2\\ & =\frac{{\mathrm{\Delta }}_i}{{\delta }_i}-c_i{\delta }_i-\frac{c_{i+1}-c_i}{3{\delta }_i}{\delta }_i^2\\ & =\frac{{\mathrm{\Delta }}_i}{{\delta }_i}-\frac{{\delta }_i}{3}(2c_i+c_{i+1})\end{array}$$

for $i=1,...,n-1$.

Substituting $d_i=\frac{c_{i+1}-c_i}{3{\delta }_i}$ and $b_i=\frac{{\mathrm{\Delta }}_i}{{\delta }_i}-\frac{{\delta }_i}{3}(2c_i+c_{i+1})$ into property 2 we can solve

$${\delta }_i{c}_i+2({\delta }_i+{\delta }_{i+1})c_{i+1}+{\delta }_{i+1}{c}_{i+2}=3(\frac{{\mathrm{\Delta }}_{i+1}}{{\delta }_{i+1}}-\frac{{\mathrm{\Delta }}_i}{{\delta }_i})$$

for $i=1,...,n-2$.

---

Property 4 Natural spline $S_1^{″}(t_1)=0$ and $S_{n-1}^{″}(t_n)=0$. (2 equations)

$$\begin{array}{rl}0& =S_1^{″}(t_1)=2c_1\\ 0& =S_{n-1}^{″}(t_n)=2c_n\end{array}$$

This gives a total of $n$ equations in $n$ unknowns $c_i$, which can be written in the matrix form

$$\left[\begin{array}{cccccc}1& 0& 0& & & \\ {\delta }_1& 2({\delta }_1+{\delta }_2)& {\delta }_2& & & \\ 0& {\delta }_2& 2({\delta }_2+{\delta }_3)& {\delta }_3& & \\ & & \ddots & \ddots & \ddots & \\ & & & {\delta }_{n-2}& 2({\delta }_{n-2}+{\delta }_{n-1})& {\delta }_{n-1}\\ & & & 0& 0& 1\end{array}\right]\left[\begin{array}{c}{c}_1\\ ⋮\\ c_n\end{array}\right]=\left[\begin{array}{c}0\\ 3(\frac{{\mathrm{\Delta }}_2}{{\delta }_2}-\frac{{\mathrm{\Delta }}_1}{{\delta }_1})\\ ⋮\\ 3(\frac{{\mathrm{\Delta }}_{n-1}}{{\delta }_{n-1}}-\frac{{\mathrm{\Delta }}_{n-2}}{{\delta }_{n-2}})\\ 0\end{array}\right].$$

---

EXAMPLE

Find the natural cubic spline through $(0,3),(1,-2)$ and $(2,1)$.

SOLUTION

${\delta }_1=1-0=1;{\delta }_2=2-1=1;{\mathrm{\Delta }}_1=-2-3=-5;{\mathrm{\Delta }}_2=1-(-2)=3$

$$\left[\begin{array}{ccc}1& 0& 0\\ 1& 4& 1\\ 0& 0& 1\end{array}\right]\left[\begin{array}{c}{c}_1\\ c_2\\ c_3\end{array}\right]=\left[\begin{array}{c}0\\ 24\\ 0\end{array}\right]$$

$[c_1,c_2,c_3]=[0,6,0]$.

$$d_1=\frac{c_2-c_1}{3{\delta }_1}=\frac{6}{3}=2;d_2=\frac{c_3-c_2}{3{\delta }_2}=\frac{-6}{3}=-2$$$$b_1=\frac{{\mathrm{\Delta }}_1}{{\delta }_1}-\frac{{\delta }_1}{3}(2c_1+c_2)=-5-\frac{1}{3}(6)=-7;b_2=\frac{{\mathrm{\Delta }}_2}{{\delta }_2}-\frac{{\delta }_2}{3}(2c_2+c_3)=3-\frac{1}{3}(12)=-1$$

Therefore, the cubic spline is

$$\begin{array}{rl}{S}_1(x)& =3-7x+0x^2+2x^3\text{ on }[0,1]\\ S_2(x)& =-2-1(x-1)+6(x-1{)}^2-2(x-1{)}^3\text{ on }[1,2].\end{array}$$