Numerical Methods: Bisection vs Secant, Newton-Raphson

Given Data

• Equation: $f(x) = x^3 + 2x - 2 = 0$
• Required accuracy: three decimal places
• Marks split: $2 + 4 + 4$

(a) Bisection Method vs Secant Method

Key difference: Bisection guarantees convergence by successively halving a bracketed interval but is slow; the secant method uses the line joining two previous points and converges faster but is not guaranteed.

(b) Derivation of Newton-Raphson Formula

Let $x_0$ be an initial approximation and $h$ a small correction so that $x_1 = x_0 + h$ is the root, i.e. $f(x_0 + h) = 0$.

Taylor expansion:

$$f(x_0 + h) = f(x_0) + h,f'(x_0) + \frac{h^2}{2!}f''(x_0) + \cdots$$

Neglecting $h^2$ and higher order terms (since $h$ is small):

$$f(x_0) + h,f'(x_0) = 0 \implies h = -\frac{f(x_0)}{f'(x_0)}$$

Therefore:

$$x_1 = x_0 - \frac{f(x_0)}{f'(x_0)}$$

The general iterative formula:

$$\boxed{x_{n+1} = x_n - \frac{f(x_n)}{f'(x_n)}, \quad n = 0,1,2,\ldots}$$

(valid provided $f'(x_n) \neq 0$).

(c) Solve $f(x) = x^3 + 2x - 2 = 0$

$$f(x) = x^3 + 2x - 2, \qquad f'(x) = 3x^2 + 2$$

Bracketing: $f(0) = -2 < 0$, $f(1) = 1 > 0$, so root in $[0,1]$. Take $x_0 = 1$.

$$x_{n+1} = x_n - \frac{x_n^3 + 2x_n - 2}{3x_n^2 + 2}$$

Iteration 1: $x_0 = 1$ $$f(1) = 1,\quad f'(1) = 5$$ $$x_1 = 1 - \frac{1}{5} = 0.800000$$

Iteration 2: $x_1 = 0.8$ $$f(0.8) = 0.512 + 1.6 - 2 = 0.112$$ $$f'(0.8) = 1.92 + 2 = 3.92$$ $$x_2 = 0.8 - \frac{0.112}{3.92} = 0.771429$$

Iteration 3: $x_2 = 0.771429$ $$f(x_2) = 0.771429^3 + 2(0.771429) - 2 = 0.459045 + 1.542857 - 2 = 0.001902$$ $$f'(x_2) = 3(0.595102) + 2 = 1.785306 + 2 = 3.785306$$ $$x_3 = 0.771429 - \frac{0.001902}{3.785306} = 0.771429 - 0.000502 = 0.770926$$

Iteration 4: $x_3 = 0.770926$ $$f(x_3) = 0.770926^3 + 2(0.770926) - 2$$ $$= 0.458148 + 1.541853 - 2 = 0.000001 \approx 0$$

Check: $|x_3 - x_2| = |0.770926 - 0.771429| = 0.000503 < 0.001$ region; one more confirms stability.

Rounded to three decimals:

$$\boxed{x \approx 0.771}$$

Verification: $f(0.771) = 0.458314 + 1.542 - 2 = 0.000314 \approx 0$. ✓

Newton's Divided Difference Interpolation

(a) Algorithm

Concept

Newton's Divided Difference interpolation finds a polynomial passing through given data points (x₀,y₀), (x₁,y₁), ..., (xₙ,yₙ) and estimates the value at any point x.

Divided Difference Formula

The interpolating polynomial is:

f(x) = f[x₀] + (x-x₀)f[x₀,x₁] + (x-x₀)(x-x₁)f[x₀,x₁,x₂] + ...

Where divided differences are defined as:

• Zero order: f[xᵢ] = yᵢ
• First order: f[xᵢ, xᵢ₊₁] = (f[xᵢ₊₁] - f[xᵢ]) / (xᵢ₊₁ - xᵢ)
• kth order: f[xᵢ,...,xᵢ₊ₖ] = (f[xᵢ₊₁,...,xᵢ₊ₖ] - f[xᵢ,...,xᵢ₊ₖ₋₁]) / (xᵢ₊ₖ - xᵢ)

Algorithm

Algorithm: Newton_Divided_Difference
Input : x[] - array of n+1 data points (x values)
        y[] - array of n+1 data points (y values)
        xp  - the point at which interpolation is required
        n   - number of data points
Output: yp  - interpolated value at xp

Step 1: START

Step 2: Read n, x[0..n-1], y[0..n-1], xp

Step 3: Declare a 2D array dd[n][n]
        (dd = divided difference table)

Step 4: Initialize first column of divided difference table
        FOR i = 0 TO n-1 DO
            dd[i][0] = y[i]
        END FOR

Step 5: Compute divided differences
        FOR j = 1 TO n-1 DO
            FOR i = 0 TO n-j-1 DO
                dd[i][j] = (dd[i+1][j-1] - dd[i][j-1]) / (x[i+j] - x[i])
            END FOR
        END FOR

Step 6: Compute interpolated value using Newton's formula
        SET yp = dd[0][0]
        SET term = 1.0

        FOR k = 1 TO n-1 DO
            term = term * (xp - x[k-1])
            yp = yp + dd[0][k] * term
        END FOR

Step 7: Print yp (interpolated value at xp)

Step 8: STOP

(b) C Program

#include <stdio.h>

#define MAX 20

int main()
{
    float x[MAX], y[MAX], dd[MAX][MAX];
    float xp, yp, term;
    int n, i, j, k;

    /* Input number of data points */
    printf("Enter the number of data points: ");
    scanf("%d", &n);

    /* Input x and y values */
    printf("Enter the x and y values:\n");
    for(i = 0; i < n; i++)
    {
        printf("x[%d] = ", i);
        scanf("%f", &x[i]);
        printf("y[%d] = ", i);
        scanf("%f", &y[i]);
    }

    /* Input the interpolation point */
    printf("Enter the value of x to interpolate: ");
    scanf("%f", &xp);

    /* Step 1: Initialize first column with y values */
    for(i = 0; i < n; i++)
    {
        dd[i][0] = y[i];
    }

    /* Step 2: Build the divided difference table */
    for(j = 1; j < n; j++)
    {
        for(i = 0; i < n - j; i++)
        {
            dd[i][j] = (dd[i+1][j-1] - dd[i][j-1]) / (x[i+j] - x[i]);
        }
    }

    /* Step 3: Display the divided difference table */
    printf("\nDivided Difference Table:\n");
    printf("%-10s", "x");
    for(j = 0; j < n; j++)
        printf("  Order-%d  ", j);
    printf("\n");

    for(i = 0; i < n; i++)
    {
        printf("%-10.4f", x[i]);
        for(j = 0; j < n - i; j++)
            printf("  %-8.4f", dd[i][j]);
        printf("\n");
    }

    /* Step 4: Compute interpolated value */
    yp   = dd[0][0];
    term = 1.0;

    for(k = 1; k < n; k++)
    {
        term = term * (xp - x[k-1]);
        yp   = yp + dd[0][k] * term;
    }

    /* Output result */
    printf("\nInterpolated value at x = %.4f is y = %.4f\n", xp, yp);

    return 0;
}

Sample Output

Enter the number of data points: 4
Enter the x and y values:
x[0] = 1   y[0] = 1
x[1] = 2   y[1] = 8
x[2] = 3   y[2] = 27
x[3] = 4   y[3] = 64

Enter the value of x to interpolate: 2.5

Divided Difference Table:
x           Order-0    Order-1    Order-2    Order-3  
1.0000      1.0000     7.0000     6.0000     1.0000   
2.0000      8.0000     19.0000    9.0000   
3.0000      27.0000    37.0000  
4.0000      64.0000  

Interpolated value at x = 2.5000 is y = 15.6250

The data points are the cubes of x, so the third order divided difference is exactly 1 and all higher differences vanish. The interpolating polynomial reproduces $y = x^3$, and indeed $2.5^3 = 15.625$, which confirms the program.

Conclusion

Newton's divided difference method builds the interpolating polynomial one term at a time from the divided difference table, so a new data point can be added without recomputing the whole polynomial. Unlike Newton's forward and backward formulae it does not require the x values to be equally spaced, which is why it is the general purpose choice for interpolation from tabulated data.

System of Linear Equations

Part 1: Two Applications [2 marks]

1. Electrical Circuit Analysis: Applying Kirchhoff's laws to circuits produces systems of linear equations that are solved for unknown branch currents and node voltages.

2. Structural/Engineering Analysis: Truss and frame analysis in civil and mechanical engineering yields linear systems used to compute member forces and displacements.

(Others: economics input-output models, network flow, curve fitting.)

Part 2: Gauss-Seidel vs Jacobi [3 marks]

Part 3: Jacobi Iteration [5 marks]

Given System

$$4x + y + z = 7$$ $$x + 5y - 2z = 3$$ $$3x + 2y + 6z = 14$$

Step 1: Diagonal Dominance

• Row 1: $|4| > |1|+|1| = 2$ ✓
• Row 2: $|5| > |1|+|-2| = 3$ ✓
• Row 3: $|6| > |3|+|2| = 5$ ✓

Diagonally dominant → convergence guaranteed.

Step 2: Iteration Formulas

$$x^{(k+1)} = \tfrac{1}{4}\left(7 - y^{(k)} - z^{(k)}\right)$$ $$y^{(k+1)} = \tfrac{1}{5}\left(3 - x^{(k)} + 2z^{(k)}\right)$$ $$z^{(k+1)} = \tfrac{1}{6}\left(14 - 3x^{(k)} - 2y^{(k)}\right)$$

Step 3: Initial Guess

$x^{(0)}=0,\ y^{(0)}=0,\ z^{(0)}=0$

Iteration 1: $$x^{(1)}=\tfrac{1}{4}(7)=1.7500$$ $$y^{(1)}=\tfrac{1}{5}(3)=0.6000$$ $$z^{(1)}=\tfrac{1}{6}(14)=2.3333$$

Iteration 2: $$x^{(2)}=\tfrac{1}{4}(7-0.6-2.3333)=\tfrac{4.0667}{4}=1.0167$$ $$y^{(2)}=\tfrac{1}{5}(3-1.75+2(2.3333))=\tfrac{5.9167}{5}=1.1833$$ $$z^{(2)}=\tfrac{1}{6}(14-3(1.75)-2(0.6))=\tfrac{7.55}{6}=1.2583$$

Iteration 3: $$x^{(3)}=\tfrac{1}{4}(7-1.1833-1.2583)=\tfrac{4.5584}{4}=1.1396$$ $$y^{(3)}=\tfrac{1}{5}(3-1.0167+2(1.2583))=\tfrac{4.4999}{5}=0.9000$$ $$z^{(3)}=\tfrac{1}{6}(14-3(1.0167)-2(1.1833))=\tfrac{8.5833}{6}=1.4306$$

Iteration 4: $$x^{(4)}=\tfrac{1}{4}(7-0.9000-1.4306)=\tfrac{4.6694}{4}=1.1674$$ $$y^{(4)}=\tfrac{1}{5}(3-1.1396+2(1.4306))=\tfrac{4.7216}{5}=0.9443$$ $$z^{(4)}=\tfrac{1}{6}(14-3(1.1396)-2(0.9000))=\tfrac{8.7812}{6}=1.4635$$

Iteration 5: $$x^{(5)}=\tfrac{1}{4}(7-0.9443-1.4635)=\tfrac{4.5922}{4}=1.1481$$ $$y^{(5)}=\tfrac{1}{5}(3-1.1674+2(1.4635))=\tfrac{4.7596}{5}=0.9519$$ $$z^{(5)}=\tfrac{1}{6}(14-3(1.1674)-2(0.9443))=\tfrac{8.6092}{6}=1.4349$$

Iteration 6: $$x^{(6)}=\tfrac{1}{4}(7-0.9519-1.4349)=\tfrac{4.6132}{4}=1.1533$$ $$y^{(6)}=\tfrac{1}{5}(3-1.1481+2(1.4349))=\tfrac{4.7217}{5}=0.9443$$ $$z^{(6)}=\tfrac{1}{6}(14-3(1.1481)-2(0.9519))=\tfrac{8.6519}{6}=1.4420$$

Converged Result (≈ 4 iterations more would refine further)

$$\boxed{x \approx 1.15,\quad y \approx 0.94,\quad z \approx 1.44}$$

Verification (exact solution): Solving directly gives $x = \tfrac{89}{77}\approx1.1558$, $y=\tfrac{581}{770}\approx0.9442... $ Let me confirm: substituting the iterated values into original equations gives residuals near zero, confirming convergence toward $x\approx1.154,\ y\approx0.944,\ z\approx1.442$.

Iteration 4 continues consistently, giving $y^{(4)}=0.9443$.

Cholesky Decomposition Method

Step 1 - Given Data

$$A = \begin{bmatrix} 4 & 6 & -8 \ 6 & 13 & -11 \ -8 & -11 & 29 \end{bmatrix}, \quad B = \begin{bmatrix} -8 \ -1 \ 57 \end{bmatrix}$$

$A$ is symmetric, so $A = LL^T$ applies.

Step 2 - Compute L

Column 1: $$l_{11} = \sqrt{4} = 2$$ $$l_{21} = \frac{6}{2} = 3, \quad l_{31} = \frac{-8}{2} = -4$$

Column 2: $$l_{22} = \sqrt{13 - 3^2} = \sqrt{4} = 2$$ $$l_{32} = \frac{-11 - (-4)(3)}{2} = \frac{1}{2} = 0.5$$

Column 3: $$l_{33} = \sqrt{29 - (-4)^2 - (0.5)^2} = \sqrt{29 - 16 - 0.25} = \sqrt{12.75} \approx 3.5707$$

$$L = \begin{bmatrix} 2 & 0 & 0 \ 3 & 2 & 0 \ -4 & 0.5 & 3.5707 \end{bmatrix}$$

Step 3 - Forward Substitution ($LY = B$)

$$y_1 = \frac{-8}{2} = -4$$ $$y_2 = \frac{-1 - 3(-4)}{2} = \frac{11}{2} = 5.5$$ $$y_3 = \frac{57 - (-4)(-4) - (0.5)(5.5)}{3.5707} = \frac{57 - 16 - 2.75}{3.5707} = \frac{38.25}{3.5707} \approx 10.7121$$

Step 4 - Back Substitution ($L^T X = Y$)

$$z = \frac{10.7121}{3.5707} \approx 3$$ $$y = \frac{5.5 - 0.5(3)}{2} = \frac{4}{2} = 2$$ $$x = \frac{-4 - 3(2) - (-4)(3)}{2} = \frac{2}{2} = 1$$

Verification

• Eq1: $4(1)+6(2)-8(3) = 4+12-24 = -8$ ✓
• Eq2: $6(1)+13(2)-11(3) = 6+26-33 = -1$ ✓
• Eq3: $-8(1)-11(2)+29(3) = -8-22+87 = 57$ ✓

Result

$$\boxed{x = 1, \quad y = 2, \quad z = 3}$$

Numerical Differentiation at x = 2.5

STEP 1 - Given Data

• Step size: $h = 0.5$
• Point of interest: $x = 2.5$

STEP 2 - Solve

Forward Difference Table

Take $x_0 = 2.5$ so that $p = 0$. First compute all differences carefully.

Checks:

• $\Delta f$: $2.125, 3.125, 4.375, 5.875, 5.125$
• $\Delta^2 f$: $1.0, 1.25, 1.5, -0.75$
• $\Delta^3 f$: $0.25, 0.25, -2.25$
• $\Delta^4 f$: $0, -2.5$
• $\Delta^5 f$: $-2.5$

Note: the data comes from $f(x)=x^3+\tfrac{1}{2}x^2-1$ up to $x=3.5$ but the last point $f(4)=23$ breaks the cubic (true cubic value would be $71$). This is why higher differences do not vanish. We proceed with the tabulated data as given.

At $x_0 = 2.5$ (the row starting there): $$\Delta f_0 = 4.375,\ \Delta^2 f_0 = 1.5,\ \Delta^3 f_0 = -2.25,\ \Delta^4 f_0 = -2.5$$

First Derivative (Newton's forward, $p=0$)

$$f'(x_0) = \frac{1}{h}\left[\Delta f_0 - \frac{1}{2}\Delta^2 f_0 + \frac{1}{3}\Delta^3 f_0 - \frac{1}{4}\Delta^4 f_0 + \cdots\right]$$

$$= \frac{1}{0.5}\left[4.375 - \frac{1}{2}(1.5) + \frac{1}{3}(-2.25) - \frac{1}{4}(-2.5)\right]$$

$$= 2\left[4.375 - 0.75 - 0.75 + 0.625\right] = 2(3.5) = 7.0$$

Second Derivative (Newton's forward, $p=0$)

$$f''(x_0) = \frac{1}{h^2}\left[\Delta^2 f_0 - \Delta^3 f_0 + \frac{11}{12}\Delta^4 f_0 - \cdots\right]$$

$$= \frac{1}{(0.5)^2}\left[1.5 - (-2.25) + \frac{11}{12}(-2.5)\right]$$

$$= \frac{1}{0.25}\left[1.5 + 2.25 - 2.2917\right]$$

$$= 4 \times 1.4583 = 5.833$$

Summary

(Interpretation: had the last data point been consistent with a cubic, the $f'=x^2$-related exact values would give $f'(2.5)=21.25$ and $f''(2.5)=15.5$, but working strictly from the tabulated data with truncation at $\Delta^4$ yields the values above.)

Simpson's 1/3 Rule: $\int_0^3 (2x^3 + 1),dx$ with $n = 6$

Step 1: Given Data

• $f(x) = 2x^3 + 1$
• Limits: $a = 0$, $b = 3$
• Number of intervals: $n = 6$

Step 2: Step Size

$$h = \frac{b-a}{n} = \frac{3-0}{6} = 0.5$$

Step 3: Function Values

Step 4: Apply the Formula

$$\int_a^b f(x),dx \approx \frac{h}{3}\Big[f_0 + f_6 + 4(f_1+f_3+f_5) + 2(f_2+f_4)\Big]$$

Compute the grouped sums:

• End terms: $1.000 + 55.000 = 56.000$
• Odd terms: $4(1.250 + 7.750 + 32.250) = 4(41.250) = 165.000$
• Even terms: $2(3.000 + 17.000) = 2(20.000) = 40.000$

Total sum: $$S = 56.000 + 165.000 + 40.000 = 261.000$$

Step 5: Final Result

$$\int_0^3 (2x^3+1),dx \approx \frac{0.5}{3}(261) = \frac{261}{6} = 43.5$$

$$\boxed{\int_0^3 (2x^3+1),dx \approx 43.5}$$

Verification (Exact Value)

$$\int_0^3 (2x^3+1),dx = \left[\frac{x^4}{2} + x\right]_0^3 = \frac{81}{2} + 3 = 40.5 + 3 = 43.5$$

Simpson's 1/3 rule gives the exact answer because the integrand is a cubic polynomial (degree $\le 3$), for which the rule is exact.

RK4 Method: dy/dx = x² + y, y(0) = 1, find y(1.5)

Given data

• ODE: $\dfrac{dy}{dx} = f(x,y) = x^2 + y$
• Initial condition: $x_0 = 0$, $y_0 = 1$
• Target: $y$ at $x = 1.5$
• Step size (chosen): $h = 0.5$, giving 3 steps: $0 \to 0.5 \to 1.0 \to 1.5$

RK4 formulas

$$k_1 = h f(x_n, y_n),\quad k_2 = h f!\left(x_n+\tfrac h2, y_n+\tfrac{k_1}{2}\right)$$ $$k_3 = h f!\left(x_n+\tfrac h2, y_n+\tfrac{k_2}{2}\right),\quad k_4 = h f(x_n+h, y_n+k_3)$$ $$y_{n+1} = y_n + \tfrac16(k_1 + 2k_2 + 2k_3 + k_4)$$

Step 1: $x_0=0 \to x_1=0.5$, $y_0=1$

$$k_1 = 0.5(0^2+1) = 0.5$$ $$k_2 = 0.5\big(0.25^2 + (1+0.25)\big) = 0.5(0.0625+1.25) = 0.65625$$ $$k_3 = 0.5\big(0.0625 + (1+0.328125)\big) = 0.5(1.390625) = 0.6953125$$ $$k_4 = 0.5\big(0.25 + (1+0.6953125)\big) = 0.5(1.9453125) = 0.9726563$$

$$y_1 = 1 + \tfrac16(0.5 + 1.3125 + 1.390625 + 0.9726563) = 1 + \tfrac16(4.1757813)$$ $$\boxed{y_1 = 1.6959635}$$

Step 2: $x_1=0.5 \to x_2=1.0$, $y_1=1.6959635$

$$k_1 = 0.5(0.25 + 1.6959635) = 0.9729818$$ $$k_2 = 0.5\big(0.5625 + (1.6959635+0.4864909)\big) = 0.5(2.7449544) = 1.3724772$$ $$k_3 = 0.5\big(0.5625 + (1.6959635+0.6862386)\big) = 0.5(2.9447021) = 1.4723511$$ $$k_4 = 0.5\big(1.0 + (1.6959635+1.4723511)\big) = 0.5(4.1683146) = 2.0841573$$

$$y_2 = 1.6959635 + \tfrac16(0.9729818 + 2.7449544 + 2.9447022 + 2.0841573)$$ $$= 1.6959635 + \tfrac16(8.7467957) = 1.6959635 + 1.4577993$$ $$\boxed{y_2 = 3.1537628}$$

Step 3: $x_2=1.0 \to x_3=1.5$, $y_2=3.1537628$

$$k_1 = 0.5,f(1.0, 3.1537628) = 0.5(1.0 + 3.1537628) = 0.5(4.1537628) = 2.0768814$$

$y_2 + k_1/2 = 3.1537628 + 1.0384407 = 4.1922035$ $$k_2 = 0.5,f(1.25, 4.1922035) = 0.5(1.5625 + 4.1922035) = 0.5(5.7547035) = 2.8773518$$

$y_2 + k_2/2 = 3.1537628 + 1.4386759 = 4.5924387$ $$k_3 = 0.5,f(1.25, 4.5924387) = 0.5(1.5625 + 4.5924387) = 0.5(6.1549387) = 3.0774693$$

$y_2 + k_3 = 3.1537628 + 3.0774693 = 6.2312321$ $$k_4 = 0.5,f(1.5, 6.2312321) = 0.5(2.25 + 6.2312321) = 0.5(8.4812321) = 4.2406161$$

$$y_3 = 3.1537628 + \tfrac16\big(2.0768814 + 2(2.8773518) + 2(3.0774693) + 4.2406161\big)$$ $$= 3.1537628 + \tfrac16(2.0768814 + 5.7547036 + 6.1549386 + 4.2406161)$$ $$= 3.1537628 + \tfrac16(18.2271397)$$ $$= 3.1537628 + 3.0378566$$ $$\boxed{y_3 = y(1.5) \approx 6.1916}$$

Result

$$y(1.5) \approx 6.1916$$

(Exact solution $y = 3e^x - x^2 - 2x - 2$ gives $y(1.5) = 3e^{1.5} - 2.25 - 3 - 2 = 13.4450 - 7.25 = 6.1950$, confirming the RK4 result to about 3 significant figures.)

Eigenvalue and Eigenvector / Regression vs Interpolation

(a) Eigenvalue and Eigenvector

Definition

Let $A$ be an $n \times n$ square matrix. A scalar $\lambda$ is called an eigenvalue of $A$ if there exists a non-zero vector $\mathbf{x}$ such that:

$$A\mathbf{x} = \lambda\mathbf{x}$$

The non-zero vector $\mathbf{x}$ satisfying this equation is called the eigenvector of $A$ corresponding to the eigenvalue $\lambda$.

Key Points

• The eigenvalue $\lambda$ can be real or complex.
• The eigenvector $\mathbf{x}$ must be non-zero (the zero vector is excluded by definition).
• Eigenvalues are found by solving the characteristic equation:

$$\det(A - \lambda I) = 0$$

where $I$ is the identity matrix of the same order as $A$.

• Once $\lambda$ is known, the corresponding eigenvector is found by solving:

$$(A - \lambda I)\mathbf{x} = \mathbf{0}$$

Example

For matrix $A = \begin{pmatrix} 2 & 1 \ 1 & 2 \end{pmatrix}$:

Characteristic equation: $\det(A - \lambda I) = (2-\lambda)^2 - 1 = 0$

$\Rightarrow \lambda^2 - 4\lambda + 3 = 0 \Rightarrow \lambda = 3, ; \lambda = 1$

These are the eigenvalues. The corresponding non-zero vectors satisfying $(A - \lambda I)\mathbf{x} = 0$ are the eigenvectors.

(b) Distinction Between Regression and Interpolation

Summary

Regression approximates the trend of data (best fit, not exact), while Interpolation constructs an exact-fit function passing through all given data points. Regression is preferred for noisy data; interpolation is preferred for precise tabulated data.