MTH603 — Final Term Summary (Lectures 23–44)

📘 Lecture 23 — Numerical Differentiation Using Difference Operators

📖 Overview: This lecture covers the numerical differentiation of functions using finite difference operators. It begins by deriving formulae for the first and second derivatives using forward, backward, and central difference operators. The lecture then introduces methods based on interpolation, Richardson's extrapolation, and ends with a transition to numerical integration, including Newton-Cotes formulae, trapezoidal rule, and Simpson's rules.

🗂️ Topics Covered

The lecture covers numerical differentiation using forward, backward, and central difference operators, including the derivation of formulae for the first and second derivatives. It proceeds to differentiation using interpolation with Newton's divided differences. Richardson's extrapolation is introduced to improve accuracy. The second half of the lecture transitions to numerical integration, covering Newton-Cotes integration formulae, the composite trapezoidal rule, and Simpson's 1/3 and 3/8 rules.

📝 Lecture Summary

Introduction

The lecture begins by reviewing the concept of an integral and its basic definition as a sum. It introduces the objective of evaluating a definite integral I = ∫f(x)dx from a to b without an analytical calculation, and discusses when this is necessary, such as when the integrand is too complicated or defined only by tabulated data. Key concepts are that integration is a summing process, and the integration methods are classified as closed forms (including endpoints) or open forms (excluding endpoints). Numerical integration approximates the integral as I = ΣΔx f(xᵢ) + Eₜ.

Newton-Cotes Integration Formulae

This method approximates the tabulated function y(x) with an interpolating polynomial Pₙ(x) of degree n and then integrates that polynomial. Using Lagrange approximation for f(x) = ΣLₖ(x)y(xₖ) with error E(x) = ∏(x) y⁽ⁿ⁺¹⁾(ξ) / (n+1)!, the integration formula becomes ∫f(x)dx ≈ Σcₖy(xₖ). The Cotes numbers (cₖ) are weighting coefficients defined by cₖ = ∫Lₖ(x)dx.

For a specific case n=1, the formula gives the Trapezoidal rule for the interval [x₀, x₁]: ∫f(x)dx = (h/2)(y₀ + y₁) - (h³/12) y''(ξ). For n=2, the integration formula yields Simpson's 1/3 rule: ∫f(x)dx = (h/3)(y₀ + 4y₁ + y₂) - (h⁵/90) y⁽ⁱᵛ⁾(ξ). This rule is based on fitting a quadratic through three points. For n=3, the formula yields Simpson's 3/8 rule: ∫f(x)dx = (3h/8)(y₀ + 3y₁ + 3y₂ + y₃) - (3h⁵/80) y⁽ⁱᵛ⁾(ξ), based on fitting a cubic through four points.

The Trapezoidal Rule (Composite Form)

To evaluate the definite integral I = ∫f(x)dx from a to b, the interval [a, b] is divided into n sub-intervals, each of size h = (b-a)/n. The composite trapezoidal rule is then derived as: ∫f(x)dx = (h/2)[y₀ + 2(y₁ + y₂ + ... + yₙ₋₁) + yₙ] + Eₙ The global error term Eₙ is given by Eₙ = - [(xₙ - x₀) h² / 12] y''(ξ). This shows the error is of the order O(h²).

Simpson's Rules (Composite Forms)

To get a composite Simpson's 1/3 rule, the interval of integration [a, b] is divided into an even number of sub-intervals, say 2N, each of width (b-a)/2N. The composite formula is: ∫f(x)dx = (h/3)[y₀ + 4(y₁ + y₃ + ... + y₂ₙ₋₁) + 2(y₂ + y₄ + ... + y₂ₙ₋₂) + y₂ₙ] + Error The global error is E = -[(x₂ₙ - x₀) h⁴ / 180] y⁽ⁱᵛ⁾(ξ), which is of order O(h⁴). The error in Simpson's 3/8 rule is of the same order, but the 1/3 rule is generally superior in magnitude.

Richardson's Extrapolation Method

This method is introduced as a way to improve the accuracy of derivative computations. Starting with a two-point formula for the derivative, the error is a series in even powers of h (E = c₁h² + c₂h⁴ + ...). By combining computed values of y'(x) using step sizes h and h/2, the h² term can be eliminated, resulting in a formula with O(h⁴) accuracy. The general formula for the extrapolation process is: Fₘ(h/2ᵐ) = [4ᵐ Fₘ₋₁(h/2ᵐ) - Fₘ₋₁(h/2ᵐ⁻¹)] / (4ᵐ - 1)

Double Integration

To evaluate a double integral of the form I = ∬f(x,y) dy dx over a rectangular region, either the trapezoidal rule or Simpson's rule is applied repeatedly. Keeping one variable fixed (e.g., y constant) and applying a rule with respect to the other variable (x) generates a table of "partial" integrals. These results are then used as integrand values for a second rule applied to the other variable (y).

⭐ Key Takeaways

• The Trapezoidal rule approximates the integral by a series of linear segments, with a global error of O(h²).
• Simpson's 1/3 rule uses a quadratic fit and has a global error of O(h⁴), making it significantly more accurate than the trapezoidal rule for the same step size. It requires an even number of intervals.
• Newton-Cotes formulae provide the theoretical framework for these integration rules by integrating an interpolating polynomial.
• The composite forms of both the Trapezoidal and Simpson's rules are practical for evaluating integrals over the entire interval [a, b] by summing the results from each sub-interval.
• Richardson's extrapolation is a powerful technique to improve the accuracy of a numerical method (like differentiation) by combining results from different step sizes, effectively increasing the order of accuracy of the method.

🧠 Quick Revision Questions

1. Derive the composite Trapezoidal rule for approximating the integral of a function f(x) from a to b and state the order of its global error.
2. What is the condition on the number of sub-intervals required for applying Simpson's 1/3 rule, and what is the formula for the composite rule?
3. How does the error magnitude of Simpson's 1/3 rule compare to the Trapezoidal rule? Explain why one is preferred.
4. Explain the basic principle behind the Newton-Cotes integration formulae.
5. How does Richardson's extrapolation improve the accuracy of a numerical derivative, and what is the order of accuracy of the result after one application?

📘 Lecture 41 — Examples of Differential Equations

📖 Overview: This lecture demonstrates the application of Euler’s Method and the Fourth-Order Runge-Kutta (RK4) Method for solving first-order differential equations numerically. It provides detailed step-by-step examples, including full calculations and algorithm implementations, showing how these methods approximate solutions to initial value problems (IVPs).

🗂️ Topics Covered

The lecture begins with a review of Euler’s Method, its recurrence relation, and algorithm. It then presents a complete numerical example using Euler’s Method to approximate the solution of a specific IVP. The lecture proceeds to recall the Fourth-Order Runge-Kutta method, provides a detailed example with all intermediate calculations for solving a differential equation from t=0 to t=0.4, and concludes with the RK4 algorithm and a computational example.

📝 Lecture Summary

Recall EULER METHOD

We considered the differential equation of first order with the initial condition y(t0) = y0. The equation is: dy/dt = f(t, y)

We obtained the solution in the form of a recurrence relation: yₘ₊₁ = yₘ + h f(tₘ, yₘ)

Euler’s method constructs wᵢ ≈ y(tᵢ) for each i = 0, 1, …, N-1 by deleting the remainder term. Thus the Euler’s Method is: w₀ = α wᵢ₊₁ = wᵢ + h f(tᵢ, wᵢ) for each i = 0, 1, …, N-1

Euler’s algorithm

The algorithm approximates the solution of the given IVP at (N+1) equally spaced numbers in the interval [a, b]. The IVP is: y′ = f(t, y), a ≤ t ≤ b, y(a) = α

INPUT: endpoints a, b; integer N; initial condition (alpha) OUTPUT: approximate w to y at the (N+1) values of t

Step 1: Set h = (b - a) / N; t = a; w = (alpha); OUTPUT (t, w) Step 2: For i = 0, 1, …, N do Step 3, 4. Step 3: Set w = w + h f(t, w); (compute wᵢ). t = a + i h (compute tᵢ) Step 4: OUTPUT (t, w) Step 5: STOP

Example (Euler’s Method)

Use Euler’s method to approximate the solution of IVP: y′ = y – t² + 1, 0 < t < 2, y(0) = 0.5 with N = 10.

Solution: Here, h = 0.2, tᵢ = 0.2i, w₀ = 0.5 and: wᵢ₊₁ = wᵢ + h (wᵢ – tᵢ² + 1) = wᵢ + 0.2[wᵢ - 0.04i² + 1] = 1.2wᵢ - 0.008i² + 0.2 for i = 0, 1, …, 9.

The exact solution is: y(t) = (t+1)² - 0.5 eᵗ

A computational example is shown using a program alg051(). For the function y-t^2+1 with endpoints 0 and 2, initial condition 0.5, and N=10, the output at t=1.000 is w=2.4581760. A second example uses the function y-3*t^2+4 with endpoints 0 and 1, initial condition 0.5, and N=10, showing output at t=1.000 as w=6.6312636.

Recall Runge-Kutta (Order Four) METHOD

The fourth-order R-K method was described as: yₙ₊₁ = yₙ + (1/6)(k₁ + 2k₂ + 2k₃ + k₄)

where: k₁ = h f(tₙ, yₙ) k₂ = h f(tₙ + h/2, yₙ + k₁/2) k₃ = h f(tₙ + h/2, yₙ + k₂/2) k₄ = h f(tₙ + h, yₙ + k₃)

Example (RK4 Method)

Solve the following differential equation: dy/dt = t + y with the initial condition y(0) = 1, using fourth-order Runge-Kutta method from t = 0 to t = 0.4 taking h = 0.1.

Solution: The fourth-order Runge-Kutta method is described as: yₙ₊₁ = yₙ + (1/6)(k₁ + 2k₂ + 2k₃ + k₄) ……….. (1)

where: k₁ = h f(tₙ, yₙ) k₂ = h f(tₙ + h/2, yₙ + k₁/2) k₃ = h f(tₙ + h/2, yₙ + k₂/2) k₄ = h f(tₙ + h, yₙ + k₃)

In this problem: f(t, y) = t + y, h = 0.1, t₀ = 0, y₀ = 1.

🔑 Definition — Step 1 (First iteration, i=0): k₁ = h f(t₀, y₀) = 0.1(1) = 0.1 k₂ = h f(t₀ + 0.05, y₀ + 0.05) = h f(0.05, 1.05) = 0.1[0.05 + 1.05] = 0.11 k₃ = h f(t₀ + 0.05, y₀ + 0.055) = 0.1(0.05 + 1.055) = 0.1105 k₄ = 0.1(0.1 + 1.1105) = 0.12105

Now, we compute from equation (1): y₁ = y₀ + (1/6)(k₁ + 2k₂ + 2k₃ + k₄) = 1 + (1/6)(0.1 + 0.22 + 0.2210 + 0.12105) = 1 + (1/6)(0.66205) = 1.11034

Therefore y(0.1) = y₁ = 1.1103

🔑 Definition — Step 2 (Second iteration, i=1): We have to find y₂ = y(0.2). We compute: k₁ = h f(t₁, y₁) = 0.1(0.1 + 1.11034) = 0.121034 k₂ = h f(t₁ + h/2, y₁ + k₁/2) = 0.1[0.15 + (1.11034 + 0.060517)] = 0.13208 k₃ = h f(t₁ + h/2, y₁ + k₂/2) = 0.1[0.15 + (1.11034 + 0.06604)] = 0.132638 k₄ = h f(t₁ + h, y₁ + k₃) = 0.1[0.2 + (1.11034 + 0.132638)] = 0.1442978

From Equation (1): y₂ = 1.11034 + (1/6)[0.121034 + 2(0.13208) + 2(0.132638) + 0.1442978] = 1.11034 + (1/6)(0.7950698) = 1.2428

🔑 Definition — Step 3 (Third iteration, i=2): We compute: k₁ = h f(t₂, y₂) = 0.1[0.2 + 1.2428] = 0.14428 k₂ = h f(t₂ + h/2, y₂ + k₁/2) = 0.1[0.25 + (1.2428 + 0.07214)] = 0.156494 k₃ = h f(t₂ + h/2, y₂ + k₂/2) = 0.1[0.3 + (1.2428 + 0.078247)] = 0.1571047 k₄ = h f(t₂ + h, y₂ + k₃) = 0.1[0.3 + (1.2428 + 0.1571047)] = 0.16999047

Using equation (1): y(0.3) = y₃ = y₂ + (1/6)(k₁ + 2k₂ + 2k₃ + k₄) = 1.399711

🔑 Definition — Step 4 (Fourth iteration, i=3): We calculate: k₁ = h f(t₃, y₃) = 0.1[0.3 + 1.3997] = 0.16997 k₂ = h f(t₃ + h/2, y₃ + k₁/2) = 0.1[0.35 + (1.3997 + 0.084985)] = 0.1834685 k₃ = h f(t₃ + h/2, y₃ + k₂/2) = 0.1[0.35 + (1.3997 + 0.091734)] = 0.1841434 k₄ = h f(t₃ + h, y₃ + k₃) = 0.1[0.4 + (1.3997 + 0.1841434)] = 0.19838434

Using them in equation (1): y(0.4) = y₄ = y₃ + (1/6)(k₁ + 2k₂ + 2k₃ + k₄) = 1.58363 which is the required result.

Runge-Kutta Order Four

w₀ = α k₁ = h f(tᵢ, wᵢ) k₂ = h f(tᵢ + h/2, wᵢ + k₁/2) k₃ = h f(tᵢ + h/2, wᵢ + k₂/2) k₄ = h f(tᵢ + h, wᵢ + k₃) wᵢ₊₁ = wᵢ + (1/6)(k₁ + 2k₂ + 2k₃ + k₄) ……….. (1)

RK4 algorithm

The algorithm approximates the solution of the given IVP at (N+1) equally spaced numbers in the interval [a, b]: y′ = f(t, y), a ≤ t ≤ b, y(a) = α

INPUT: endpoints a, b; integer N; initial condition (alpha) OUTPUT: approximate w to y at the (N+1) values of t

Step 1: Set h = (b - a) / N; t = a; w = (alpha); OUTPUT (t, w) Step 2: For i = 0, 1, …, N do Step 3 - 5. Step 3: K₁ = h f(t, w) K₂ = h f(t + h/2, w + K₁/2) K₃ = h f(t + h/2, w + K₂/2) K₄ = h f(t + h, w + K₃) w = w + (1/6)(k₁ + 2k₂ + 2k₃ + k₄)

Set w = w + h f(t, w); (compute wᵢ). t = a + i h (compute tᵢ) Step 4: OUTPUT (t, w) Step 5: STOP

A computational example is shown using a program alg052(). For the function y-t^2+1 with endpoints 0 and 2, initial condition 0.5, and N=10, the output at t=2.000 is w=5.3053630.

⭐ Key Takeaways

The most critical points from this lecture are: Euler’s Method and the Fourth-Order Runge-Kutta Method are numerical techniques for approximating solutions to first-order differential equations. Euler’s Method uses a simple recurrence relation wᵢ₊₁ = wᵢ + h f(tᵢ, wᵢ) and is less accurate. The Fourth-Order Runge-Kutta Method uses a weighted average of four slope estimates (k₁, k₂, k₃, k₄) and is significantly more accurate. Both methods require the user to input the function f(t,y), endpoints, initial condition, and number of subintervals N. The step size h is computed as (b-a)/N, and the algorithms generate approximate values wᵢ at equally spaced points tᵢ.

🧠 Quick Revision Questions

1. What is the recurrence relation for Euler’s Method for solving dy/dt = f(t, y) with initial condition y(t₀) = y₀?
2. In the Fourth-Order Runge-Kutta method, what are the four slope estimates (k₁, k₂, k₃, k₄) and how are they combined to compute yₙ₊₁?
3. In the Euler’s Method example with y′ = y – t² + 1, y(0) = 0.5, N=10, what is the step size h and what is the resulting recurrence relation for wᵢ₊₁?
4. In the RK4 example with dy/dt = t + y, y(0) = 1, h = 0.1, what is the value of y(0.1) after the first iteration?
5. What are the INPUT and OUTPUT specifications for both the Euler’s Method algorithm and the RK4 algorithm?

📘 Lecture 42 — Examples of Numerical Differentiation

📖 Overview: This lecture covers practical applications of numerical differentiation, including the simplest differentiation formula, three-point and five-point formulas, and their error analysis. It then transitions into numerical integration methods, covering the Trapezoidal rule, Simpson's rule, closed and open Newton-Cotes formulas, composite numerical integration, improper integrals, and an introduction to using MAPLE for numerical computation.

🗂️ Topics Covered

The lecture begins with examples of numerical differentiation using the two-point forward difference formula, three-point formulas (endpoint and midpoint), and a five-point formula, all with error analysis and worked examples. It then covers numerical integration with the Trapezoidal and Simpson's rules, closed and open Newton-Cotes formulas, composite numerical integration for improved accuracy, improper integrals handled via Taylor series expansion and variable substitution, and concludes with a comprehensive introduction to MAPLE software, its syntax, and its capabilities for symbolic and numerical computation.

📝 Lecture Summary

Examples of Numerical Differentiation

The simplest formula for differentiation is given by f'(x₀) = [f(x₀ + h) - f(x₀)]/h - (h/2)f''(ξ), where ξ is between x₀ and x₀+h.

🔑 Definition — Forward Difference Formula: f'(x₀) ≈ [f(x₀+h) - f(x₀)]/h with error term -(h/2)f''(ξ).

📌 Example: Let f(x)= ln x and x₀=1.8. Using the forward difference quotient [f(1.8+h)-f(1.8)]/h to approximate f'(1.8) = 1/1.8 ≈ 0.555556. The error bound is |h·f''(ξ)/2| ≤ h/(2·(1.8)²). For h=0.1, approximation=0.5406722; h=0.01, approximation=0.5540180; h=0.001, approximation=0.5554013. The exact value is 0.555556.

When nodes are equally spaced with x₁=x₀+h and x₂=x₀+2h, the three-point formulas become especially useful.

🔑 Definition — Three-Point Endpoint Formula: f'(x₀) = [ -3f(x₀) + 4f(x₀+h) - f(x₀+2h) ]/(2h) + (h²/3)f⁽³⁾(ξ₀), where ξ₀ lies between x₀ and x₀+2h.

🔑 Definition — Three-Point Midpoint Formula: f'(x₀) = [ f(x₀+h) - f(x₀-h) ]/(2h) - (h²/6)f⁽³⁾(ξ₁), where ξ₁ lies between x₀-h and x₀+h.

📌 Example: Given table values for f(x)=xeˣ at x=1.8,1.9,2.0,2.1,2.2. Exact f'(2.0) = (x+1)eˣ|₂.₀ = 22.167168. Using three-point formulas: For h=0.1 (endpoint): [ -3(14.778112)+4(17.148957)-19.855030 ]/0.2 = 22.032310. For h=-0.1 (endpoint): [ -3(14.778112)+4(12.703199)-10.889365 ]/(-0.2) = 22.0054525. For h=0.1 (midpoint): [ 17.148957-12.703199 ]/0.2 = 22.228790. For h=0.2 (midpoint): [ 19.855030-10.889365 ]/0.4 = 22.414163.

🔑 Definition — Five-Point Midpoint Formula: f'(x₀) = [ f(x₀-2h) - 8f(x₀-h) + 8f(x₀+h) - f(x₀+2h) ]/(12h).

📌 Example: Using five-point formula with h=0.1: f'(2.0) = [ f(1.8) - 8f(1.9) + 8f(2.1) - f(2.2) ]/1.2 = [10.889365 - 8(12.703199) + 8(17.148957) - 19.855030]/1.2 = 22.166999. Errors: endpoint formulas ≈ 1.35×10⁻¹ and 1.13×10⁻¹, midpoint formulas ≈ -6.16×10⁻² and -2.47×10⁻¹, five-point formula ≈ 1.69×10⁻⁴. Clearly, the five-point formula gives the superior result.

📌 Example: Approximating f'(0.900) for f(x)=sin x (true value cos(0.900)=0.62161) using [ f(0.900+h)-f(0.900-h) ]/(2h). For h=0.001: approx=0.62500, error=0.00339; h=0.002: approx=0.62250, error=0.00089; h=0.005: approx=0.62200, error=0.00039; h=0.010: approx=0.62150, error=-0.00011; h=0.020: approx=0.62150, error=-0.00011; h=0.050: approx=0.62140, error=-0.00021; h=0.100: approx=0.62055, error=-0.00106.

Examples of Numerical Integration

The Trapezoidal rule for a function f on [0,2] is ∫₀² f(x)dx ≈ (h/2)[f(x₀)+f(x₁)] = f(0)+f(2). Simpson's rule for f on [0,2] is ∫₀² f(x)dx ≈ (h/3)[f(x₀)+4f(x₁)+f(x₂)] = (1/3)[f(0)+4f(1)+f(2)].

📌 Example: Table of exact vs. approximate integrals for various functions on [0,2]:

Common closed Newton-Cotes formulas with error terms:

• n=1: Trapezoidal rule: ∫x₀ˣ¹ f(x)dx = (h/2)[f(x₀)+f(x₁)] - (h³/12)f''(ξ), x₀<ξ<x₁
• n=2: Simpson's rule: ∫x₀ˣ² f(x)dx = (h/3)[f(x₀)+4f(x₁)+f(x₂)] - (h⁵/90)f⁽⁴⁾(ξ), x₀<ξ<x₂
• n=3: ∫x₀ˣ³ f(x)dx = (3h/8)[f(x₀)+3f(x₁)+3f(x₂)+f(x₃)] - (3h⁵/80)f⁽⁴⁾(ξ)
• n=4: ∫x₀ˣ⁴ f(x)dx = (2h/45)[7f(x₀)+32f(x₁)+12f(x₂)+32f(x₃)+7f(x₄)] - (8h⁷/945)f⁽⁶⁾(ξ)

🔑 Definition — Midpoint Rule (n=0 open formula): ∫x₋₁ˣ¹ f(x)dx = 2hf(x₀) + (h³/3)f''(ξ), where x₋₁<ξ<x₁.

📌 Example: Approximating ∫₀^(π/4) sin x dx = 1 - √2/2:

Composite Numerical Integration

📌 Example 1: Approximating ∫₀^π sin x dx with absolute error < 0.00002 using Composite Simpson's rule. Error bound: πh⁴/180 = π⁵/(180n⁴) < 0.00002, giving n ≥ 18. Using n=20: h=π/20, ∫₀^π sin x dx ≈ (π/60)[2∑ⱼ₌₁⁹ sin(jπ/10) + 4∑ⱼ₌₁¹⁰ sin((2j-1)π/20)] = 2.000006. Exact value is 2. Using Composite Trapezoidal rule requires n ≥ 360 for same accuracy. With n=20, Composite Trapezoidal gives 1.9958860, clearly insufficient.

📌 Industrial Application: Probability that a toilet roll has ≥250 sheets: P(y≥250) = ∫₂₅₀²⁷⁰ 0.3515e^(-0.3881(y-252.2)²) dy. a) Using single segment Trapezoidal rule: I ≈ (b-a)[f(a)+f(b)]/2 = 20(0.053721+1.3888×10⁻⁵⁴)/2 = 0.53721 b) True error (using Maple exact value 0.97377): Eₜ = 0.97377 - 0.53721 = 0.43656 c) Absolute relative true error: |εₜ| = |0.97377-0.53721|/0.97377 × 100 = 44.832%

Improper Integrals

📌 Example: Approximating ∫₀¹ eˣ/√x dx using Composite Simpson's rule with h=0.25. Use fourth Taylor polynomial P₄(x)=1+x+x²/2+x³/6+x⁴/24. Then ∫₀¹ P₄(x)/√x dx = lim_{M→0⁺} [2√x + (2/3)x^(3/2) + (1/5)x^(5/2) + (1/21)x^(7/2) + (1/108)x^(9/2)]_M¹ ≈ 2.9235450. Define G(x)=[eˣ-P₄(x)]/√x for 0<x≤1 and G(0)=0. Apply Composite Simpson's: (0.25/3)[0+4(0.0000170)+2(0.0004013)+4(0.0026026)+0.0099485]=0.0017691. Thus ∫₀¹ eˣ/√x dx ≈ 2.9235450+0.0017691=2.9253141, with error bounded by 0.0000217.

📌 Example: Approximating I = ∫₁^∞ x^(-3/2) sin(1/x) dx. Make substitution t=x⁻¹ to get I = ∫₀¹ t^(1/2) sin t dt. Use P₄(t)=t - t³/6. Then I = ∫₀¹ [sin t - t + t³/6]/t^(1/2) dt + ∫₀¹ [t^(1/2) - t^(7/2)/6] dt = ∫₀¹ [sin t - t + t³/6]/t^(1/2) dt + [(2/3)t^(3/2) - (1/21)t^(9/2)]₀¹ = integral term + 0.61904761. Apply Composite Simpson's with n=8 to remaining integral: I = 0.0014890097 + 0.61904761 = 0.62053661, accurate to within 4.0×10⁻⁸.

Introduction to MAPLE

Maple is a comprehensive computer system for advanced mathematics including interactive algebra, calculus, discrete mathematics, graphics, and numerical computation. It provides a unique environment for rapid development of mathematical programs.

Key syntax elements: ; processes and shows output; : processes and hides output; := assigns values; +, -, *, / arithmetic; ^ for power; sqrt() for square root; evalf() for floating-point evaluation; I for imaginary unit; Pi for π; % for last output; %% for second-to-last output.

Maple is case-sensitive. x*y gives product, xy is one variable. Use exp(1) for constant e.

Built-in data capabilities include arbitrary-precision floating-point numbers controlled by the Digits variable. Default Digits=10. Use evalf(value, N) for temporary precision.

💡 Why this matters: Maple handles both symbolic and numerical integration. Use int(f, x) for symbolic; evalf(Int(f, x=a..b)) for numerical without symbolic first. The inert Int prevents symbolic evaluation.

Basic plotting: plot(x^2, x=-2..2), plot([x, x^2, x^3], x=-2..2), plot3d(4-x^2-y^2, x=-3..3, y=-2..2).

Eigenvalues and eigenvectors: Eigenvals(A, vecs) or eigenvectors(A) from linalg package.

Maple can solve: non-linear equations (Newton's, Bisection), systems of linear equations, numerical integration, ODEs.

Numerical integration functions: evalf(Int(f, x=a..b)), evalf(Int(f, a..b)), alg041() for Simpson's method, alg044() for double integrals with Simpson's, alg045() for Gaussian quadrature for double integrals.

⭐ Key Takeaways

This lecture demonstrates that numerical differentiation accuracy improves significantly from forward difference (error O(h)) to three-point formulas (O(h²)) to five-point formulas (O(h⁴)), as shown by errors of 10⁻¹ to 10⁻⁴ for the xeˣ example. For numerical integration, Simpson's rule with error O(h⁵) is vastly superior to the Trapezoidal rule with O(h³), requiring n=18 vs n=360 for the same accuracy when integrating sin x. Composite formulas dramatically improve accuracy over single-segment formulas. Improper integrals can be handled by subtracting Taylor series expansion to remove singularities, then applying standard quadrature. MAPLE is a powerful computational environment that handles both symbolic and numerical mathematics, with special attention needed for its syntax and the distinction between symbolic (int) and numerical (evalf(Int)) integration.

🧠 Quick Revision Questions

1. What is the error term for the three-point midpoint formula for numerical differentiation?
2. Why does the five-point formula give better results than the three-point formulas in the xeˣ example?
3. What is the error term for Simpson's rule (n=2 closed Newton-Cotes formula)?
4. How many subintervals n are required for Composite Simpson's rule to approximate ∫₀^π sin x dx with error < 0.00002?
5. In MAPLE, what is the difference between int(f, x=a..b) and evalf(Int(f, x=a..b))?

📘 Lecture 44 — Solution of Non-Linear Equations

📖 Overview: This lecture presents numerical methods for solving non-linear equations and systems of linear equations. It covers root-finding techniques including Bisection, Regula-Falsi, Newton-Raphson, and Muller's methods, alongside direct and iterative methods for linear systems, interpolation formulas, numerical differentiation, integration, and solutions to differential equations.

🗂️ Topics Covered

The lecture covers methods for solving non-linear equations (Bisection, Regula-Falsi, Method of Iteration, Newton-Raphson, Muller's, Graeffe's Root Squaring), systems of linear equations (Gaussian Elimination, Gauss-Jordan, Crout's Reduction, Jacobi, Gauss-Seidel, Relaxation), eigenvalue problems (Power Method), interpolation and finite difference operators, numerical differentiation and integration (Trapezoidal Rule, Double Integration), and differential equations (Taylor Series, Euler, Runge-Kutta, Predictor-Corrector methods).

📝 Lecture Summary

Newton-Raphson Method

An approximation to the root is given by the formula ( x = x_0 - \frac{f(x_0)}{f'(x_0)} ). Better and successive approximations ( x_2, x_3, ..., x_n ) to the root are obtained from the N-R Formula: ( x_{n+1} = x_n - \frac{f(x_n)}{f'(x_n)} ).

Newton's algorithm to find a solution to ( f(x) = 0 ):

1. Input initial approximation ( p_0 ); tolerance TOL; maximum number of iterations ( N_0 ).
2. Set ( i = 1 ).
3. While ( i < N_0 ) do Steps 4-7.
4. Set ( p = p_0 - f(p_0) / f'(p_0) ) (compute ( p_i )).
5. If ( \text{Abs}(p - p_0) < \text{TOL} ), OUTPUT ( p ); STOP.
6. Set ( i = i + 1 ).
7. Set ( p_0 = p ) (Update ( p_0 )).
8. OUTPUT "The method failed after ( N_0 ) iterations". STOP.

Example: Using Maple to solve ( \cos(x) - x = 0 ). The Maple command fsolve(cos(x)-x) yields an approximate solution of 0.73908513. With initial approximation 0.7853981635, tolerance 0.00005, and 25 max iterations, Newton's method took 3 iterations, with F(P) ≈ 0.0000000000.

Another Example: ( \sin(x) - 1 = 0 ). With initial approximation 0.17853, the method took 15 iterations to converge to 1.57076292 with F(P) = 6.0000000000e-10.

🔑 Definition — Newton-Raphson Method: An iterative method for finding successively better approximations to the roots of a real-valued function ( f(x) = 0 ), using the formula ( x_{n+1} = x_n - f(x_n)/f'(x_n) ).

📐 Formula: ( x_{n+1} = x_n - \frac{f(x_n)}{f'(x_n)} ) → The new approximation is found by subtracting the ratio of the function value to its derivative at the current point.

📌 Example: For ( \cos(x) - x = 0 ), starting from ( x_0 = 0.7853981635 ), successive iterations gave ( x_1 = 0.739536134 ), ( x_2 = 0.739085178 ), and ( x_3 = 0.739085133 ), achieving tolerance in 3 iterations.

💡 Why this matters: Newton-Raphson converges quadratically (very fast) when the initial guess is close to the true root, making it one of the most widely used root-finding algorithms in numerical analysis.

Bisection Method

Example: Using Maple's alg021() to solve ( x^3 + 4x^2 - 10 = 0 ). Input endpoints A=1, B=2, tolerance 0.0005, max 25 iterations. After 11 iterations, the approximate solution was 1.36474609 with F(P) = -0.00798926.

Another Bisection Example: ( \cos(x) = 0 ). Input endpoints 1 and 2. After 11 iterations, the approximate solution was 1.57080078 with F(P) = -0.00000445.

🔑 Definition — Bisection Method (Bolzano): A root-finding method that repeatedly bisects an interval and selects a subinterval in which a root must lie, based on the Intermediate Value Theorem.

Method of False Position (Regula-Falsi)

The first approximation to the root of ( f(x) = 0 ) is given by ( x_{n+1} = x_n - \frac{x_n - x_{n-1}}{f(x_n) - f(x_{n-1})} f(x_n) ). Here ( f(x_{n-1}) ) and ( f(x_{n+1}) ) are of opposite sign.

Example: Using alg025() to solve ( \cos(x) - x = 0 ). Input endpoints ( P_0 = 0.5, P_1 = 0.7853981635 ). After 4 iterations, approximate solution was 0.73908486 with F(P) = 0.00000045.

🔑 Definition — Regula-Falsi Method: A root-finding algorithm that uses a linear interpolation between two points with opposite function signs to approximate the root, converging faster than bisection in many cases.

📐 Formula: ( x_{n+1} = x_n - \frac{x_n - x_{n-1}}{f(x_n) - f(x_{n-1})} f(x_n) ) → The next approximation is found by linear interpolation between the last two points.

Method of Iteration

This method can be applied by rewriting ( f(x) = 0 ) in the form ( x = \phi(x) ). Successive approximations are: ( x_1 = \phi(x_0), x_2 = \phi(x_1), ..., x_{n+1} = \phi(x_n) ).

N-R Formula and Secant Method

The Secant method formula is ( x_{n+1} = \frac{x_{n-1}f(x_n) - x_n f(x_{n-1})}{f(x_n) - f(x_{n-1})} ), which converges faster than linear but slower than Newton's quadratic convergence.

Muller's Method

A better approximation is obtained by ( x_{i+1} = x_i + h_i \lambda ), where ( \lambda = \frac{-2 f_i \delta_i}{g \pm [g^2 - 4 f_i \delta_i \lambda (f_{i-2} \lambda - f_{i-1} \delta_i + f_i)]^{1/2}} ), with ( \lambda = h_i / h_{i-1} ), ( h_i = x_i - x_{i-1} ), and ( \delta_i = 1 + \lambda_i ).

Systems of Linear Equations

Gaussian Elimination

The solution is obtained in two stages:

1. The given system is reduced to an equivalent upper triangular form using elementary transformations.
2. The upper triangular system is solved using back substitution.

Example: Using alg061() for a 4-equation system. The reduced system output showed a triangular matrix, and the solution vector was [-7.00000000, 3.00000000, 2.00000000, 2.00000000] with 1 row interchange.

Gauss-Jordan Elimination

A variation of Gaussian method where elements above and below the diagonal are simultaneously made zero.

Crout's Reduction Method

The coefficient matrix [A] is decomposed into [L][U], where [L] is a lower-triangular matrix and [U] is an upper-triangular matrix with 1's on its main diagonal. Example: for a 3×3 matrix, ( [L][U] = [A] ) where L has non-zero elements on and below the diagonal, and U has 1's on the diagonal with non-zero elements above.

Jacobi's Iterative Method

An iterative method where the ((r+1))th approximation is given by: ( x_i^{(r+1)} = \frac{b_i}{a_{ii}} - \sum_{j \neq i} \frac{a_{ij}}{a_{ii}} x_j^{(r)} )

No element of ( x_i^{(r+1)} ) replaces ( x_i^{(r)} ) entirely for the next cycle.

Example: Using alg071(), after 10 iterations with tolerance 0.001, the solution vector was 1.00011860, 1.99976795, -0.99982814, 0.99978598.

Gauss-Seidel Method

Corresponding elements of ( x_i^{(r+1)} ) replace those of ( x_i^{(r)} ) as soon as they become available, also called the method of Successive Displacement.

Relaxation Method

Another iterative method due to Southwell, where ( dx_i = R_i / a_{ii} ).

Eigenvalue Problems

Power Method

The result: ( u^{(k)} = [A] v^{(k-1)} = q_k v^{(k)} ), where ( q_k ) is the desired largest eigenvalue and ( v^{(k)} ) is the corresponding eigenvector.

Interpolation — Finite Difference Operators

Forward Difference Operator (Δ)

( \Delta y_x = y_{x+h} - y_x = f(x+h) - f(x) ) ( \Delta^2 y_x = \Delta y_{x+h} - \Delta y_x )

Backward Difference Operator (∇)

( \nabla y_x = y_x - y_{x-h} = f(x) - f(x-h) ) ( \nabla^k y_i = \nabla^{k-1} y_i - \nabla^{k-1} y_{i-1}, \quad i = n, (n-1), \ldots, k )

Central Difference Operator (δ)

( \delta y_x = y_{x+(h/2)} - y_{x-(h/2)} = f(x + h/2) - f(x - h/2) ) ( \delta^k y_i = \delta^{k-1} y_{i+1/2} - \delta^{k-1} y_{i-1/2} )

Shift Operator (E)

( E f(x) = f(x + h) ), ( E^n f(x) = f(x + nh) ), ( E^{-1} f(x) = f(x - h) ), ( E^{-n} f(x) = f(x - nh) )

Average Operator (μ)

( \mu f(x) = \frac{1}{2} [f(x + h/2) + f(x - h/2)] = \frac{1}{2} [y_{x+(h/2)} + y_{x-(h/2)}] )

Differential Operator (D)

( D f(x) = \frac{d}{dx} f(x) = f'(x) ), ( D^2 f(x) = f''(x) )

Important Relations

( \Delta = E - 1 ), ( \nabla = 1 - E^{-1} = \frac{E-1}{E} ), ( \delta = E^{1/2} - E^{-1/2} ), ( hD = \log E ), ( \mu = \frac{1}{2} (E^{1/2} + E^{-1/2}) )

Newton's Forward Difference Interpolation Formula

( f(x_0 + ph) = f(x_0) + p\Delta f(x_0) + \frac{p(p-1)}{2!} \Delta^2 f(x_0) + \frac{p(p-1)(p-2)}{3!} \Delta^3 f(x_0) + \cdots + \frac{p(p-1)\cdots(p-n+1)}{n!} \Delta^n f(x_0) + \text{Error} )

Newton's Backward Difference Interpolation Formula

( f(x_n + ph) = f(x_n) + p\nabla f(x_n) + \frac{p(p+1)}{2!} \nabla^2 f(x_n) + \frac{p(p+1)(p+2)}{3!} \nabla^3 f(x_n) + \cdots + \frac{p(p+1)(p+2)\cdots(p+n-1)}{n!} \nabla^n f(x_n) + \text{Error} )

Here ( p = \frac{x - x_n}{h} ).

Lagrange's Interpolation Formula

( y = f(x) = \sum_{i=0}^n \left( \prod_{j \neq i} \frac{x - x_j}{x_i - x_j} \right) y_i )

Newton's Divided Difference Interpolation Formula

( y = f(x) = y_0 + (x - x_0) y[x_0, x_1] + (x - x_0)(x - x_1) y[x_0, x_1, x_2] + \cdots + (x - x_0)(x - x_1)\cdots(x - x_{n-1}) y[x_0, x_1, \ldots, x_n] )

The first order divided difference: ( y[x_0, x_1] = \frac{y_1 - y_0}{x_1 - x_0} )

Numerical Differentiation and Integration

D expressed in terms of Δ: ( D = \frac{1}{h} \left( \Delta - \frac{\Delta^2}{2} + \frac{\Delta^3}{3} - \frac{\Delta^4}{4} + \frac{\Delta^5}{5} - \cdots \right) )

Using backward difference: ( hD = -\log(1 - \nabla) ), so ( D = \frac{1}{h} \left( \nabla + \frac{\nabla^2}{2} + \frac{\nabla^3}{3} + \frac{\nabla^4}{4} + \cdots \right) )

Using central difference: ( D = \frac{1}{h} \left( \delta - \frac{1}{24} \delta^3 + \frac{3}{640} \delta^5 - \cdots \right) )

Richardson's Extrapolation: ( F_m \left( \frac{h}{2^m} \right) = \frac{4^m F_{m-1}(h/2^m) - F_{m-1}(h/2^{m-1})}{4^m - 1} ), for ( m = 1, 2, 3, \ldots )

Trapezoidal Rule

( \int_{x_0}^{x_n} f(x) dx = \frac{h}{2} (y_0 + 2y_1 + 2y_2 + \cdots + 2y_{n-1} + y_n) + E_n )

The integral ( \int_{x_0}^{x_3} f(x) dx = \frac{3h}{8} (y_0 + 3y_1 + 3y_2 + y_3) - \frac{3h^5}{80} y^{(iv)}(\xi) )

Double Integration

The integral ( I = \int \left[ \int (x, y) dx \right] dy ) represents volume under a surface.

Differential Equations

Taylor Series Method

( y(t) = y(t_0) + (t - t_0) y'(t_0) + \frac{(t - t_0)^2}{2!} y''(t_0) + \frac{(t - t_0)^3}{3!} y'''(t_0) + \frac{(t - t_0)^4}{4!} y^{IV}(t_0) + \cdots )

Euler Method

Solution in the form of a recurrence relation: ( y_{m+1} = y_m + h f(t_m, y_m) )

Modified Euler's Method: ( y_{m+1} = y_m + h \left[ \frac{f(t_m, y_m) + f(t_{m+1}, y_{m+1}^{(1)})}{2} \right] )

Fourth-Order Runge-Kutta Method

( y_{n+1} = y_n + \frac{1}{6} (k_1 + 2k_2 + 2k_3 + k_4) ) where:

• ( k_1 = h f(t_n, y_n) )
• ( k_2 = h f(t_n + h/2, y_n + k_1/2) )
• ( k_3 = h f(t_n + h/2, y_n + k_2/2) )
• ( k_4 = h f(t_n + h, y_n + k_3) )

Milne's Predictor-Corrector Pair

Predictor: ( y_{n+1} = y_{n-3} + \frac{4h}{3} (2y'{n-2} - y'{n-1} + 2y'n) ) Corrector: ( y{n+1} = y_{n-1} + \frac{h}{3} (y'_{n-1} + 4y'n + y'{n+1}) )

Adam's Predictor Formula

( y_{n+1} = y_n + \frac{h}{24} (55 f_n - 59 f_{n-1} + 37 f_{n-2} - 9 f_{n-3}) + \frac{251}{720} h \nabla^4 f_n ) Alternatively: ( y_{n+1} = y_n + \frac{h}{24} [55 y'n - 59 y'{n-1} + 37 y'{n-2} - 9 y'{n-3}] + \frac{251}{720} h \nabla^4 y'_n )

⭐ Key Takeaways

1. Newton-Raphson is the fastest root-finding method (quadratic convergence) but requires the derivative and a good initial guess; the formula is ( x_{n+1} = x_n - f(x_n)/f'(x_n) ). 2. For solving linear systems, Gaussian Elimination reduces to upper triangular form followed by back substitution, while Jacobi and Gauss-Seidel are iterative methods differing in whether updates are immediate. 3. The five finite difference operators (Δ, ∇, δ, E, μ) are linked by fundamental relations like ( \Delta = E - 1 ), ( \nabla = 1 - E^{-1} ), and ( hD = \log E ), enabling numerical differentiation and interpolation. 4. Key interpolation formulas include Newton's Forward/Backward (for equally spaced data) and Lagrange's formula (for unequally spaced data), with divided differences generalizing Newton interpolation to arbitrary spacing. 5. For differential equations, the fourth-order Runge-Kutta method balances accuracy and simplicity, while Predictor-Corrector methods (Milne, Adams) offer higher accuracy through iterative refinement.

🧠 Quick Revision Questions

1. What is the Newton-Raphson iteration formula and what condition is necessary for the method to converge quickly?
2. How does the Gauss-Seidel method differ from Jacobi's method in iterative linear system solving?
3. State the relationship between the forward difference operator Δ and the shift operator E.
4. Write the Newton's forward difference interpolation formula and identify what variable p represents.
5. What are the four k-values computed in the fourth-order Runge-Kutta method, and how are they combined to produce ( y_{n+1} )?