BSc CSIT (TU) Science Numerical Method (BSc CSIT, CSC207) Question Paper 2078 Nepal

Curve Fitting

Curve fitting is the process of constructing a curve (a mathematical function) that has the best fit to a series of observed data points $(x_i, y_i)$. The goal is to find a smooth functional relationship $y = f(x)$ that approximately represents the trend of the data, so it can be used for interpolation, extrapolation, and prediction.

Least Squares Method for a Straight Line

We wish to fit the line

to $n$ data points $(x_i, y_i),\ i = 1, 2, \dots, n$.

For each point the residual (error) is

The least squares principle requires choosing $a$ and $b$ so that the sum of the squares of the residuals is minimum:

Deriving the Normal Equations

$S$ is minimum when its partial derivatives with respect to $a$ and $b$ vanish:

Simplifying gives the two normal equations:

Solving for the Coefficients

Solving these simultaneous equations:

where $\bar{x} = \frac{1}{n}\sum x_i$ and $\bar{y} = \frac{1}{n}\sum y_i$.

Procedure

1. Form a table of $x_i,\ y_i,\ x_i^2,\ x_i y_i$ and compute $\sum x_i,\ \sum y_i,\ \sum x_i^2,\ \sum x_i y_i$.
2. Substitute into the normal equations.
3. Solve for $a$ and $b$.
4. The required fitted line is $y = a + bx$.

The resulting line passes through the centroid $(\bar{x}, \bar{y})$ and gives the best straight-line approximation in the least-squares sense.

Finite Difference Method (FDM) for PDEs

The finite difference method solves a partial differential equation by replacing the continuous derivatives with difference approximations evaluated on a discrete grid (mesh) of points covering the solution domain.

Steps:

1. Cover the region with a rectangular grid of mesh size $h$ in $x$ and $k$ in $y$, giving nodes $(x_i, y_j)$.
2. Replace each partial derivative in the PDE by a finite-difference approximation at the interior nodes.
3. This converts the PDE into a system of algebraic equations relating the nodal values $u_{i,j}$.
4. Impose the boundary conditions, then solve the system (directly or by iteration such as Gauss-Seidel / Liebmann's method).

Standard central-difference approximations (mesh size $h$):

Five-Point Formula for Laplace's Equation

Laplace's equation in two dimensions is

Using a square grid ($h = k$) and substituting the central-difference approximations above:

Multiplying through by $h^2$:

Solving for the central value gives the standard five-point (diagonal) formula:

Interpretation: the value of $u$ at any interior node equals the average of the values at its four neighbouring nodes (left, right, above, below). The system formed at all interior nodes is solved together with the prescribed boundary values, usually iteratively by Liebmann's (Gauss-Seidel) method.

Bisection Method

The bisection method is a bracketing (interval-halving) technique to find a real root of $f(x) = 0$ where $f$ is continuous.

Principle (Intermediate Value Theorem): If $f(a)\cdot f(b) < 0$, then a root lies between $a$ and $b$.

Algorithm:

1. Choose $a, b$ such that $f(a)\,f(b) < 0$.
2. Compute the midpoint $c = \dfrac{a+b}{2}$.
3. If $f(c) = 0$ (or $|f(c)|$ is within tolerance), $c$ is the root.
4. If $f(a)\,f(c) < 0$, the root lies in $[a, c]$, so set $b = c$; otherwise set $a = c$.
5. Repeat until $|b - a|$ (or $|f(c)|$) is less than the required accuracy.

The method always converges (linearly) since the interval halves each step.

Root of $x^3 - x - 1 = 0$

Let $f(x) = x^3 - x - 1$.

$f(1) = 1 - 1 - 1 = -1 < 0$ and $f(2) = 8 - 2 - 1 = 5 > 0$, so a root lies in $[1, 2]$.

The successive midpoints converge to $1.3247$.

(The exact root is $1.32472\ldots$.)

Newton's Divided Difference Interpolation

This formula constructs an interpolating polynomial through data points $(x_0, y_0), (x_1, y_1), \dots, (x_n, y_n)$ with unequally spaced $x$-values.

Divided differences are defined recursively:

• First order: $f[x_0, x_1] = \dfrac{f(x_1) - f(x_0)}{x_1 - x_0}$
• Second order: $f[x_0, x_1, x_2] = \dfrac{f[x_1, x_2] - f[x_0, x_1]}{x_2 - x_0}$
• In general: $f[x_0, \dots, x_k] = \dfrac{f[x_1, \dots, x_k] - f[x_0, \dots, x_{k-1}]}{x_k - x_0}$

Newton's divided difference formula:

Procedure: Build a divided-difference table, take the leading (top) divided differences as coefficients, and substitute the value of $x$ at which the function is required.

Advantages: Works for unequal spacing; new data points can be added without recomputing earlier terms (unlike Lagrange's method).

Power Method

The power method is an iterative technique to find the dominant (largest in magnitude) eigenvalue $\lambda_1$ and its corresponding eigenvector of a square matrix $A$.

Principle: If $A$ has a unique dominant eigenvalue $\lambda_1$ ($|\lambda_1| > |\lambda_2| \ge \cdots$), then repeatedly multiplying an arbitrary starting vector by $A$ makes it converge to the eigenvector of $\lambda_1$.

Algorithm:

1. Choose an initial nonzero vector $x_0$ (e.g. $[1, 1, \dots, 1]^T$).
2. Compute $y_{k+1} = A x_k$.
3. Let $\lambda^{(k+1)} =$ the largest-magnitude component of $y_{k+1}$.
4. Normalize: $x_{k+1} = \dfrac{1}{\lambda^{(k+1)}}\, y_{k+1}$.
5. Repeat steps 2-4 until successive $\lambda$ values (and vectors) agree to the desired accuracy.

At convergence, $\lambda^{(k)} \to \lambda_1$ (the largest eigenvalue) and $x_k$ is the corresponding eigenvector.

Notes: Convergence is linear and faster when $|\lambda_2/\lambda_1|$ is small. The inverse power method (applied to $A^{-1}$) gives the smallest eigenvalue.

Fitting a Second-Degree Parabola (Least Squares)

We fit

to $n$ data points $(x_i, y_i)$. The residual for each point is $e_i = y_i - (a + bx_i + cx_i^2)$, and we minimize

Normal Equations

Setting $\dfrac{\partial S}{\partial a} = \dfrac{\partial S}{\partial b} = \dfrac{\partial S}{\partial c} = 0$ gives three normal equations:

Procedure

1. Tabulate $x, y, x^2, x^3, x^4, xy, x^2y$ and compute the required sums.
2. Substitute the sums into the three normal equations.
3. Solve the simultaneous linear equations (e.g. by elimination) for $a, b, c$.
4. The fitted parabola is $y = a + bx + cx^2$.

This yields the best second-degree polynomial fit in the least-squares sense.

Types of Error

Let $X$ be the true value and $X'$ the approximate (computed) value.

• Absolute error: the magnitude of the difference,

• Relative error: absolute error per unit true value,

• Percentage error: relative error expressed as a percentage,

Sources of Errors in Numerical Computation

1. Inherent (input/data) errors – errors already present in the data due to limitations of measuring instruments or imperfect physical assumptions.
2. Round-off errors – caused by representing numbers with a finite number of digits (finite word length of the computer), e.g. retaining only a few decimal places.
3. Truncation errors – caused by approximating an infinite process by a finite one, e.g. terminating a Taylor/infinite series after a few terms.
4. Modelling / formulation errors – arise when the mathematical model is a simplification of the real problem.
5. Human/blunders – mistakes in coding, data entry, or arithmetic.

Round-off and truncation errors are the two principal sources controllable in numerical methods.

Secant Method

The secant method finds a root of $f(x)=0$ using two previous approximations. The next estimate is the $x$-intercept of the secant line joining $(x_{n-1}, f(x_{n-1}))$ and $(x_n, f(x_n))$:

Algorithm:

1. Choose two initial guesses $x_0, x_1$.
2. Compute $x_2$ using the formula above.
3. Replace $(x_0, x_1) \to (x_1, x_2)$ and repeat until $|x_{n+1} - x_n|$ is within tolerance.

This is essentially the Newton-Raphson method with the derivative $f'(x_n)$ replaced by the difference quotient $\dfrac{f(x_n)-f(x_{n-1})}{x_n - x_{n-1}}$.

Comparison with Newton-Raphson

Summary: The secant method avoids computing derivatives and is cheaper per iteration, but converges a little slower than the faster (quadratic) Newton-Raphson method.

False Position (Regula Falsi) Method

The false position method is a bracketing method that, like bisection, requires two points $a, b$ with $f(a)\,f(b) < 0$ so that a root lies between them. Instead of taking the midpoint, it takes the point where the chord (straight line) joining $(a, f(a))$ and $(b, f(b))$ crosses the x-axis.

Iteration formula:

Algorithm:

1. Choose $a, b$ with $f(a)\,f(b) < 0$.
2. Compute $x$ from the formula above.
3. If $f(a)\,f(x) < 0$, set $b = x$; else set $a = x$ (keep the sub-interval that still brackets the root).
4. Repeat until $|f(x)|$ (or successive $x$ values) is within the desired accuracy.

Geometric Interpretation

The two points $(a, f(a))$ and $(b, f(b))$ lie on opposite sides of the x-axis. The straight chord connecting them intersects the x-axis at $x$. This intersection $x$ is taken as the improved approximation to the root. Replacing the curve $y=f(x)$ locally by this chord (a linear interpolation) is why it is called the method of false position: as iterations proceed, the chord's intercept approaches the actual root where the curve cuts the x-axis.

Convergence is always guaranteed (the root stays bracketed) and is usually faster than bisection, though one endpoint may remain fixed.

Lagrange's Interpolation Formula

For $n+1$ data points $(x_0, y_0), (x_1, y_1), \dots, (x_n, y_n)$ with unequally (or equally) spaced abscissae, the interpolating polynomial is

where the Lagrange basis polynomials are

For three points the formula is

Each $L_i(x) = 1$ at $x = x_i$ and $0$ at every other node, so the polynomial passes exactly through all the data points.

Example

Given $(0, 1), (1, 3), (2, 7)$, find $y$ at $x = 1.5$.

So $y(1.5) \approx 4.75$.

Trapezoidal Rule — Derivation

To evaluate $\displaystyle I = \int_{a}^{b} f(x)\,dx$, divide $[a, b]$ into $n$ equal sub-intervals of width $h = \dfrac{b-a}{n}$, with nodes $x_0=a, x_1, \dots, x_n=b$ and ordinates $y_i = f(x_i)$.

Single strip: Over each sub-interval the curve is approximated by a straight line, so the area of the strip $[x_i, x_{i+1}]$ is that of a trapezium:

Composite rule: Summing over all $n$ strips:

That is: (end ordinates) + 2 × (sum of remaining ordinates), all multiplied by $h/2$.

Error Term

The error in the composite trapezoidal rule is

Since the error is proportional to $h^2$, the trapezoidal rule is a second-order method and is exact for polynomials of degree $\le 1$ (straight lines).

Gauss-Seidel Iterative Method

The Gauss-Seidel method is an iterative technique for solving a system of $n$ linear equations $AX = B$. It refines an initial guess of the unknowns until convergence.

Consider the system:

Step 1 — Rearrange each equation to make one variable the subject (solve the $i$-th equation for $x_i$ using the diagonal element):

Step 2 — Iterate: Start with an initial guess (e.g. $x_1=x_2=x_3=0$). In each sweep, compute $x_1$, then immediately use the new $x_1$ when computing $x_2$, and use the new $x_1, x_2$ when computing $x_3$. This use of the most recently updated values distinguishes it from the Jacobi method and speeds up convergence.

Step 3 — Stop when the values stop changing within the required tolerance, i.e. $|x_i^{(k+1)} - x_i^{(k)}| < \varepsilon$ for all $i$.

Convergence condition: The method is guaranteed to converge if the coefficient matrix is diagonally dominant, i.e.

so the equations should be arranged to make the largest coefficients lie on the diagonal before iterating.