Numerical Methods

Equations: $2x+2y+z=6$, $4x+2y+3z=4$, $x-y=-1$.

Augmented matrix:

$\left[\begin{array}{ccc|c}2&2&1&6\\4&2&3&4\\1&-1&0&-1\end{array}\right]$

Pivot 1. The largest entry in column 1 is $4$, so swap $R_1\leftrightarrow R_2$:

$\left[\begin{array}{ccc|c}4&2&3&4\\2&2&1&6\\1&-1&0&-1\end{array}\right]$

$R_2\to R_2-0.5R_1$ and $R_3\to R_3-0.25R_1$:

$\left[\begin{array}{ccc|c}4&2&3&4\\0&1&-0.5&4\\0&-1.5&-0.75&-2\end{array}\right]$

Pivot 2. The larger entry in column 2 is $-1.5$, so swap $R_2\leftrightarrow R_3$:

$\left[\begin{array}{ccc|c}4&2&3&4\\0&-1.5&-0.75&-2\\0&1&-0.5&4\end{array}\right]$

$R_3\to R_3+\tfrac23R_2$:

$\left[\begin{array}{ccc|c}4&2&3&4\\0&-1.5&-0.75&-2\\0&0&-1&\tfrac83\end{array}\right]$

Back substitution

• $z=-\dfrac83$
• $-1.5y-0.75\left(-\dfrac83\right)=-2\Rightarrow -1.5y=-4\Rightarrow y=\dfrac83$
• $4x+2\left(\dfrac83\right)+3\left(-\dfrac83\right)=4\Rightarrow 4x=\dfrac{20}{3}\Rightarrow x=\dfrac53$

$x=\dfrac53\approx1.667,\; y=\dfrac83\approx2.667,\; z=-\dfrac83\approx-2.667$

Check in the third equation: $\frac53-\frac83=-1$, which is correct.

Use the five-point formula with $h=1$:

$f_{i+1,j}+f_{i-1,j}+f_{i,j+1}+f_{i,j-1}-4f_{i,j}=h^2\cdot2x^2y^2$

The interior points are $(1,1),(2,1),(1,2),(2,2)$ and all boundary values are $0$. Let $f_{1,1}=p$, $f_{2,1}=f_{1,2}=a$ and $f_{2,2}=q$ (the problem is symmetric in $x$ and $y$).

• At $(1,1)$: $2a-4p=2(1)^2(1)^2=2$
• At $(2,1)$: $p+q-4a=2(2)^2(1)^2=8$
• At $(2,2)$: $2a-4q=2(2)^2(2)^2=32$

From the first, $p=\dfrac{a-1}{2}$. From the third, $q=\dfrac{a-16}{2}$. Substituting into the second:

$\frac{2a-17}{2}-4a=8\;\Rightarrow\;-3a=16.5\;\Rightarrow\;a=-5.5$

Then $p=-3.25$ and $q=-10.75$.

$f(1,1)=-3.25,\; f(2,1)=f(1,2)=-5.5,\; f(2,2)=-10.75$.

Use the five-point finite-difference formula with $h=1$:

$f_{i+1,j}+f_{i-1,j}+f_{i,j+1}+f_{i,j-1}-4f_{i,j}=h^2\cdot2x_i^2y_j^2$

The interior points are $(1,1),(2,1),(1,2),(2,2)$, with $f=0$ on the boundary. By symmetry in $x$ and $y$, write $f_{1,1}=p$, $f_{2,1}=f_{1,2}=a$, $f_{2,2}=q$.

• Point $(1,1)$: $2a-4p=2$
• Point $(2,1)$: $p+q-4a=8$
• Point $(2,2)$: $2a-4q=32$

From the first and third equations, $p=\dfrac{a-1}{2}$ and $q=\dfrac{a-16}{2}$. Substituting into the second:

$\frac{2a-17}{2}-4a=8\;\Rightarrow\;-3a=16.5\;\Rightarrow\;a=-5.5$

So $p=-3.25$ and $q=-10.75$.

$f(1,1)=-3.25,\; f(2,1)=f(1,2)=-5.5,\; f(2,2)=-10.75$.

Augmented matrix:

$\left[\begin{array}{ccc|c} 2 & 3 & 1 & 9 \\ 1 & 4 & 5 & 15 \\ 3 & 1 & 2 & 10 \end{array}\right]$

Pivot 1. The largest entry in column 1 is 3 (row 3). Swap $R_1 \leftrightarrow R_3$:

$\left[\begin{array}{ccc|c} 3 & 1 & 2 & 10 \\ 1 & 4 & 5 & 15 \\ 2 & 3 & 1 & 9 \end{array}\right]$

Eliminate with $R_2 \leftarrow R_2 - \tfrac13 R_1$ and $R_3 \leftarrow R_3 - \tfrac23 R_1$:

$\left[\begin{array}{ccc|c} 3 & 1 & 2 & 10 \\ 0 & 11/3 & 13/3 & 35/3 \\ 0 & 7/3 & -1/3 & 7/3 \end{array}\right]$

Pivot 2. Between $11/3$ and $7/3$, the larger is $11/3$ (already in row 2), so no swap is needed. Apply $R_3 \leftarrow R_3 - \tfrac{7}{11}R_2$:

$\left[\begin{array}{ccc|c} 3 & 1 & 2 & 10 \\ 0 & 11/3 & 13/3 & 35/3 \\ 0 & 0 & -34/11 & -56/11 \end{array}\right]$

Back substitution.

$z = \frac{56/11}{34/11} = \frac{28}{17} \approx 1.6471$

$\frac{11}{3}y = \frac{35}{3} - \frac{13}{3}\cdot\frac{28}{17} \;\Rightarrow\; 11y = 35 - \frac{364}{17} = \frac{231}{17} \;\Rightarrow\; y = \frac{21}{17} \approx 1.2353$

$3x = 10 - y - 2z = 10 - \frac{21}{17} - \frac{56}{17} = \frac{93}{17} \;\Rightarrow\; x = \frac{31}{17} \approx 1.8235$

Solution: $x \approx 1.8235,\; y \approx 1.2353,\; z \approx 1.6471$.

Check in equation 1: $2(1.8235) + 3(1.2353) + 1.6471 = 3.647 + 3.706 + 1.647 = 9.000$.

Derivation of Simpson's 3/8 rule. Fit a cubic through four equally spaced points $x_0, x_1=x_0+h, x_2, x_3$ with values $y_0,\dots,y_3$. Put $x=x_0+ph$, so $dx=h\,dp$ and $p$ runs from $0$ to $3$. Newton's forward formula gives

$y(p)=y_0+p\,\Delta y_0+\frac{p(p-1)}{2}\Delta^2y_0+\frac{p(p-1)(p-2)}{6}\Delta^3y_0$

Integrating term by term over $p=0$ to $3$:

$\int_0^3 dp=3,\quad \int_0^3 p\,dp=\frac92,\quad \int_0^3\frac{p(p-1)}{2}dp=\frac94,\quad \int_0^3\frac{p(p-1)(p-2)}{6}dp=\frac38$

$I=h\left[3y_0+\frac92\Delta y_0+\frac94\Delta^2y_0+\frac38\Delta^3y_0\right]$

Substitute $\Delta y_0=y_1-y_0$, $\Delta^2y_0=y_2-2y_1+y_0$, $\Delta^3y_0=y_3-3y_2+3y_1-y_0$ and collect terms:

$I=\int_{x_0}^{x_3}f(x)\,dx\approx\frac{3h}{8}\left[y_0+3y_1+3y_2+y_3\right]$

Secant method for $f(x)=x^2-4x-10$. The iteration is

$x_{n+1}=x_n-\frac{f(x_n)(x_n-x_{n-1})}{f(x_n)-f(x_{n-1})}$

With $x_1=4$ and $x_2=2$: $f(4)=-10$ and $f(2)=-14$.

The root is $x\approx 5.7417$, which agrees with the exact value $2+\sqrt{14}=5.7417$.

Boundary value problem. A differential equation of order two or higher in which the conditions are given at more than one point, usually at the two ends of the interval (for example $y(a)=\alpha$ and $y(b)=\beta$), is a boundary value problem. In an initial value problem, all conditions are given at a single point.

Shooting method. Replace the BVP by an IVP. Guess the missing slope $y'(0)=s$, integrate to $x=b$, and compare $y(b)$ with the required value. Adjust $s$ by linear interpolation (secant) until $y(b)$ matches the boundary value.

Here $y''=6x^2$, $y(0)=1$, $y(1)=2$. Write $z=y'$, so $y'=z$ and $z'=6x^2$. Integrating the second equation gives $z=2x^3+s$ and then

$y(x)=\frac{x^4}{2}+sx+1$

This is exact, so it agrees with Euler/RK stepping on this problem.

• Trial 1, $s=0$: $y(1)=0.5+0+1=1.5$.
• Trial 2, $s=1$: $y(1)=0.5+1+1=2.5$.

The target is $y(1)=2$. Because $y(1)$ is linear in $s$,

$s=0+\frac{2-1.5}{2.5-1.5}(1-0)=0.5$

Then with $h=0.5$:

$y(0.5)=\frac{0.5^4}{2}+0.5(0.5)+1=0.03125+0.25+1=1.28125$

Answer: $y(0.5)\approx 1.28125$.

Algorithm (Lagrange interpolation)

1. Read $n$, the data points $(x_i,y_i)$ for $i=0,\dots,n-1$, and the point $x_p$ at which to interpolate.
2. Set $sum=0$.
3. For $i=0$ to $n-1$:
   • set $term=y_i$;
   • for $j=0$ to $n-1$ with $j\ne i$, set $term=term\times\dfrac{x_p-x_j}{x_i-x_j}$;
   • set $sum=sum+term$.
4. Print $sum$ as $f(x_p)$.
5. Stop.

The formula used is

$f(x_p)=\sum_{i=0}^{n-1}y_i\prod_{\substack{j=0\\ j\ne i}}^{n-1}\frac{x_p-x_j}{x_i-x_j}$

C program

#include <stdio.h>

int main(void) {
    int n, i, j;
    double x[20], y[20], xp, sum = 0.0, term;

    printf("Enter number of data points: ");
    scanf("%d", &n);
    printf("Enter x and y values:\n");
    for (i = 0; i < n; i++)
        scanf("%lf %lf", &x[i], &y[i]);

    printf("Enter the value of x to interpolate: ");
    scanf("%lf", &xp);

    for (i = 0; i < n; i++) {
        term = y[i];
        for (j = 0; j < n; j++) {
            if (j != i)
                term *= (xp - x[j]) / (x[i] - x[j]);
        }
        sum += term;
    }

    printf("f(%.4f) = %.6f\n", xp, sum);
    return 0;
}

The x values must be distinct, otherwise the program divides by zero.

Newton-Raphson iteration. For $f(x)=0$,

$x_{n+1}=x_n-\frac{f(x_n)}{f'(x_n)}$

Let $\alpha$ be a simple root, so $f(\alpha)=0$ and $f'(\alpha)\ne0$. Define the error $e_n=x_n-\alpha$.

Expand $f(\alpha)$ about $x_n$ by Taylor's theorem, with $\xi$ between $x_n$ and $\alpha$:

$0=f(\alpha)=f(x_n)+(\alpha-x_n)f'(x_n)+\frac{(\alpha-x_n)^2}{2}f''(\xi)$

Divide by $f'(x_n)$ and rearrange:

$x_n-\frac{f(x_n)}{f'(x_n)}-\alpha=\frac{f''(\xi)}{2f'(x_n)}(x_n-\alpha)^2$

The left side is $x_{n+1}-\alpha=e_{n+1}$, so

$e_{n+1}=\frac{f''(\xi)}{2f'(x_n)}\,e_n^2$

As $x_n\to\alpha$, we have $\xi\to\alpha$ and

$|e_{n+1}|\approx\left|\frac{f''(\alpha)}{2f'(\alpha)}\right|\,|e_n|^2=C\,|e_n|^2$

The new error is proportional to the square of the previous error, so the order of convergence is $2$ (quadratic). The number of correct digits roughly doubles with every iteration once $x_n$ is close to the root.

Let $X=e^{t/4}$. Then $T=bX+a$ is linear in $X$, and least squares applies.

Normal equations:

$4a+7.76803\,b=377,\qquad 7.76803\,a+16.23775\,b=775.493$

Solving gives $b\approx37.63$ and $a\approx21.17$, so

$T=37.63\,e^{t/4}+21.17$

At $t=6$, $e^{1.5}=4.4817$:

$T=37.63(4.4817)+21.17=168.64+21.17\approx189.8$

Estimated temperature at $t=6$ min is about $190\,^\circ$C.

Divided difference table

Each entry is $\dfrac{\text{difference of the two entries to its left}}{x_{i+k}-x_i}$. For example, $\dfrac{19-7}{3-1}=6$.

Newton's divided difference polynomial (the 4th difference is $0$, as expected for a cubic):

$f(x)=0+7(x-1)+6(x-1)(x-2)+1(x-1)(x-2)(x-3)$

At $x=1.5$:

$f(1.5)=7(0.5)+6(0.5)(-0.5)+(0.5)(-0.5)(-1.5)=3.5-1.5+0.375=2.375$

$f(1.5)=2.375$, which equals $1.5^3-1$ exactly.

Characteristic equation $|A-\lambda I|=0$. Expanding along the third row:

$(3-\lambda)\big[(1-\lambda)(2-\lambda)-6\big]=(3-\lambda)(\lambda^2-3\lambda-4)=(3-\lambda)(\lambda-4)(\lambda+1)=0$

Eigenvalues: $\lambda=4,\;3,\;-1$.

For $\lambda=4$: $(A-4I)v=0$ gives $-3x+6y+z=0$, $x-2y=0$, $-z=0$. So $z=0$, $x=2y$:

$v_1=\begin{bmatrix}2\\1\\0\end{bmatrix}$

For $\lambda=3$: $-2x+6y+z=0$, $x-y=0$, $0=0$. So $x=y$ and $z=-4x$:

$v_2=\begin{bmatrix}1\\1\\-4\end{bmatrix}$

For $\lambda=-1$: $2x+6y+z=0$, $x+3y=0$, $4z=0$. So $z=0$, $x=-3y$:

$v_3=\begin{bmatrix}3\\-1\\0\end{bmatrix}$

Any non-zero scalar multiple of each vector is also an eigenvector.

Step size $h=1$ s. Use central differences at $t=7$ with $s(6)=14.5$, $s(7)=19.5$, $s(8)=25.5$.

Velocity

$v=\frac{ds}{dt}\approx\frac{s(8)-s(6)}{2h}=\frac{25.5-14.5}{2}=5.5\ \text{m/s}$

Acceleration

$a=\frac{d^2s}{dt^2}\approx\frac{s(8)-2s(7)+s(6)}{h^2}=\frac{25.5-39.0+14.5}{1}=1.0\ \text{m/s}^2$

At $t=7$ s: velocity $\approx5.5$ m/s and acceleration $\approx1.0$ m/s$^2$.

$I=\int_2^4(x^4+1)dx$ with $n=8$, so $h=\dfrac{4-2}{8}=0.25$. The rule is

$I\approx\frac h2\Big[f_0+f_8+2(f_1+\dots+f_7)\Big]$

Sum of interior values $=669.2656$.

$I\approx0.125\big[17+257+2(669.2656)\big]=0.125(1612.5312)\approx201.566$

$I\approx201.57$. The exact value is $200.4$, so the trapezoidal rule overestimates slightly because $f$ is convex.

$y'=3x^2+1$, $y(1)=2$, $h=0.2$. Euler's formula: $y_{n+1}=y_n+h\,f(x_n,y_n)$.

$y(2)\approx9.12$. The exact solution is $y=x^3+x$, giving $y(2)=10$, so Euler with this step size underestimates.

True error: $E_t=\text{true value}-\text{approximate value}$.

Relative (true) error: $\varepsilon_t=\dfrac{E_t}{\text{true value}}$, often given as a percentage.

Bisection method (derivation). Suppose $f$ is continuous on $[a,b]$ and $f(a)f(b)<0$. By the intermediate value theorem there is a root in $(a,b)$. Take the midpoint

$c=\frac{a+b}{2}$

• If $f(c)=0$, then $c$ is the root.
• If $f(a)f(c)<0$, the root lies in $[a,c]$, so set $b=c$.
• Otherwise the root lies in $[c,b]$, so set $a=c$.

Repeat. The interval halves every step, so after $n$ steps its width is $(b-a)/2^n$.

Application: $f(x)=2x^3-2x-5$ on $[1,2]$. Here $f(1)=-5<0$ and $f(2)=7>0$.

After 10 iterations the root is $x\approx1.6006$ (accurate to about three decimal places).

Applications of interpolation

• Estimating values between tabulated data, for example from a table of sines, logarithms or steam tables.
• Reconstructing missing readings from experiments or sensors.
• Numerical differentiation and integration, which are built on interpolating polynomials.
• Computer graphics and image resizing (curves, surfaces, resampling).

Interpolation versus regression

Lagrange interpolation at $x=10$ with $x=5,6,9,11$ and $y=13,14,15,16$:

$L_0=\frac{(10-6)(10-9)(10-11)}{(5-6)(5-9)(5-11)}=\frac{-4}{-24}=\frac16$

$L_1=\frac{(10-5)(10-9)(10-11)}{(6-5)(6-9)(6-11)}=\frac{-5}{15}=-\frac13$

$L_2=\frac{(10-5)(10-6)(10-11)}{(9-5)(9-6)(9-11)}=\frac{-20}{-24}=\frac56$

$L_3=\frac{(10-5)(10-6)(10-9)}{(11-5)(11-6)(11-9)}=\frac{20}{60}=\frac13$

$f(10)=13\left(\tfrac16\right)+14\left(-\tfrac13\right)+15\left(\tfrac56\right)+16\left(\tfrac13\right)=\tfrac{13}{6}-\tfrac{28}{6}+\tfrac{75}{6}+\tfrac{32}{6}=\tfrac{92}{6}$

$f(10)\approx15.33$. The weights sum to $\frac16-\frac13+\frac56+\frac13=1$, as they must.

Numerical integration means approximating a definite integral $\int_a^b f(x)\,dx$ by a weighted sum of function values, $\sum w_if(x_i)$. It is needed when $f$ has no closed-form antiderivative or when $f$ is known only as a table of values.

Application: finding the distance travelled from tabulated speed readings, or the area of an irregular plot of land from measured offsets.

Multi-segment trapezoidal rule. Divide $[a,b]$ into $n$ equal segments of width $h=\dfrac{b-a}{n}$:

$I\approx\frac h2\Big[f(x_0)+2\sum_{i=1}^{n-1}f(x_i)+f(x_n)\Big]$

Algorithm

1. Read $a$, $b$ and $n$.
2. Set $h=(b-a)/n$ and $sum=f(a)+f(b)$.
3. For $i=1$ to $n-1$: set $x=a+ih$ and $sum=sum+2f(x)$.
4. Set $I=(h/2)\times sum$.
5. Print $I$ and stop.

C program (the example integrand is $f(x)=\dfrac1{1+x^2}$)

#include <stdio.h>

double f(double x) {
    return 1.0 / (1.0 + x * x);
}

int main(void) {
    double a, b, h, sum, x, I;
    int n, i;

    printf("Enter lower limit a, upper limit b: ");
    scanf("%lf %lf", &a, &b);
    printf("Enter number of segments n: ");
    scanf("%d", &n);

    h = (b - a) / n;
    sum = f(a) + f(b);
    for (i = 1; i < n; i++) {
        x = a + i * h;
        sum += 2.0 * f(x);
    }
    I = (h / 2.0) * sum;

    printf("Integral = %.6f\n", I);
    return 0;
}

Data: $x=1,2,4,8,10$ and $y=0,1,5,21,27$.

Divided difference table

For example, $f[2,4,8]=\dfrac{4-2}{8-2}=\dfrac13$ and $f[2,4,8,10]=\dfrac{-1/6-1/3}{10-2}=-\dfrac1{16}$.

Interpolating polynomial

$P(x)=(x-1)+\tfrac13(x-1)(x-2)-\tfrac1{144}(x-1)(x-2)(x-4)(x-8)$

First derivative at $x=2$. Let $g(x)=(x-1)(x-2)(x-4)(x-8)$. Since $(x-2)$ is a factor, $g'(2)=(2-1)(2-4)(2-8)=12$. Then

$P'(2)=1+\tfrac13(2\cdot2-3)-\tfrac{12}{144}=1+\tfrac13-\tfrac1{12}=\tfrac{17}{12}\approx1.417$

Second derivative at $x=2$. Write $g=(x-2)h$ with $h=(x-1)(x-4)(x-8)=x^3-13x^2+44x-32$. Then $g''(2)=2h'(2)=2(3\cdot4-26\cdot2+44)=2(4)=8$, and

$P''(2)=\tfrac23-\tfrac8{144}=\tfrac23-\tfrac1{18}=\tfrac{11}{18}\approx0.611$

$f'(2)\approx1.417$ and $f''(2)\approx0.611$.

Least squares method. Given data $(x_i,y_i)$, choose the curve's parameters so that the sum of squares of the vertical deviations $S=\sum[y_i-f(x_i)]^2$ is minimum. Setting the partial derivative of $S$ with respect to each parameter to zero gives the normal equations. For $y=a_0+a_1x+a_2x^2$:

$\sum y=na_0+a_1\sum x+a_2\sum x^2$ $\sum xy=a_0\sum x+a_1\sum x^2+a_2\sum x^3$ $\sum x^2y=a_0\sum x^2+a_1\sum x^3+a_2\sum x^4$

Sums for the data ($n=7$): $\sum x=28$, $\sum x^2=140$, $\sum x^3=784$, $\sum x^4=4676$, $\sum y=60$, $\sum xy=294$, $\sum x^2y=1634$.

$7a_0+28a_1+140a_2=60$ $28a_0+140a_1+784a_2=294$ $140a_0+784a_1+4676a_2=1634$

Solving gives $a_0=\dfrac87=1.1429$, $a_1=\dfrac{73}{42}=1.7381$ and $a_2=\dfrac1{42}=0.0238$.

$y\approx1.1429+1.7381x+0.0238x^2$