MTH603 — Midterm Summary (Lectures 1–22)

📘 Lecture 1 — Numerical Analysis

📖 Overview: This introductory lecture covers the fundamental concepts of number systems used in computing and the representation of numbers on computers. It is essential because numerical analysis relies on understanding how numbers are stored and manipulated in binary form, and introduces the course contents including solution of equations, interpolation, differentiation, integration, and differential equations.

🗂️ Topics Covered

The lecture begins with an introduction to numerical analysis and problem solving using computers, then covers number systems including decimal, binary, octal, and hexadecimal bases. It explains how real numbers are represented in different bases, provides examples of converting decimal numbers to binary (both integers and fractions), and demonstrates conversion from binary to octal. The course contents are outlined at the start.

📝 Lecture Summary

Introduction

This chapter begins with basic concepts of representation of numbers on computers and errors introduced during computation. Problem solving using computers and the steps involved are also discussed in brief.

Number (s) System (s)

In our daily life, we use numbers based on the decimal system, which uses ten symbols 0, 1,...,9 and the number 10 is called the base of the system. When a base N is given, we need N different symbols 0, 1, 2, ..., (N – 1) to represent an arbitrary number.

The number systems commonly used in computers are:

An arbitrary real number, a, can be written as: a = aₘNᵐ + aₘ₋₁Nᵐ⁻¹ + ... + a₁N¹ + a₀ + a₋₁N⁻¹ + ... + a₋ₘN⁻ᵐ

In the binary system, it has the form: a = aₘ2ᵐ + aₘ₋₁2ᵐ⁻¹ + ... + a₁2¹ + a₀ + a₋₁2⁻¹ + ... + a₋ₘ2⁻ᵐ

The decimal number 1729 is represented and calculated as: (1729)₁₀ = 1×10³ + 7×10² + 2×10¹ + 9×10⁰

The decimal equivalent of binary number 10011001 is: 1×2⁰ + 0×2⁻¹ + 0×2⁻² + 1×2⁻³ + 1×2⁻⁴ + 0×2⁻⁵ + 0×2⁻⁶ + 1×2⁻⁷ = 1 + 1/8 + 1/16 + 1/128 = (1.1953125)₁₀

Electronic computers use the binary system whose base is 2. The two symbols used are 0 and 1, which are called binary digits or simply bits. The internal representation of any data within a computer is in binary form. However, data input and output of numerical results are preferred in decimal system. Within the computer, arithmetic is carried out in binary form.

🔑 Definition — Base: The number of distinct symbols used in a number system; for decimal it is 10, for binary it is 2. 🔑 Definition — Bit: A binary digit, either 0 or 1, the fundamental unit of data in computers.

Conversion of decimal number 47 into its binary equivalent

The conversion is done by repeated division by 2:

(47)₁₀ = (101111)₂

📐 Formula: Decimal to binary integer conversion → Repeatedly divide by 2, reading remainders from last (most significant bit) to first (least significant bit). 📌 Example: (47)₁₀ → Divide by 2 repeatedly; remainders from last to first give 101111.

Binary equivalent of the decimal fraction 0.7625

The conversion is done by repeated multiplication by 2:

Reading the integer parts from top to bottom: (0.7625)₁₀ = (0.110000110...)₂

📐 Formula: Decimal to binary fraction conversion → Repeatedly multiply by 2, recording integer parts in order. 📌 Example: (0.7625)₁₀ → Multiply by 2 repeatedly; integer parts in order give 0.110000110...

Conversion (59)₁₀ into binary and then into octal

First, convert 59 to binary:

(59)₁₀ = (111011)₂

Then convert binary to octal by grouping bits in threes from right: (111011)₂ = 111 011 = (73)₈

📐 Formula: Binary to octal conversion → Group binary digits in threes from right, convert each group to its octal equivalent. 📌 Example: (111011)₂ → Group as 111 and 011 → (7)(3) = (73)₈

💡 Why this matters: Understanding number system conversions is critical for numerical analysis because all computer arithmetic is performed in binary, and errors can arise during conversion between bases.

⭐ Key Takeaways

The lecture establishes that numerical analysis involves representing and solving mathematical problems on computers, where numbers are stored in binary form. Four number systems are commonly used in computing: decimal (base 10), binary (base 2), octal (base 8), and hexadecimal (base 16). Converting integers from decimal to binary requires repeated division by 2, reading remainders from last to first. Converting decimal fractions to binary requires repeated multiplication by 2, recording integer parts in order. Binary numbers can be converted to octal by grouping bits in threes, and understanding these conversions is fundamental to managing computational errors in numerical methods.

🧠 Quick Revision Questions

1. What is the binary equivalent of the decimal number 47?
2. How do you convert a decimal fraction to binary, and what is the binary equivalent of 0.7625?
3. What is the octal representation of the binary number 111011?
4. Why is the binary system used in electronic computers instead of the decimal system?
5. What are bits and why are they important in numerical computing?

📘 Lecture 2 — Errors in Computations

📖 Overview: This lecture introduces the fundamental types of errors that arise in numerical computations. Understanding these errors is critical because numerically computed solutions are never exact, and knowing the sources and growth of errors helps in selecting appropriate algorithms and interpreting results with confidence.

🗂️ Topics Covered

The lecture classifies errors into three main categories: inherent errors, which exist in the problem statement before any computation; local round-off errors, which arise due to finite computer word length and binary-decimal conversion; and local truncation errors, which result from approximating infinite series by a finite number of terms. Each error type is defined, explained with examples, and illustrated with practical significance.

📝 Lecture Summary

Errors in Computations

Numerically computed solutions are subject to certain errors. It may be fruitful to identify the error sources and their growth while classifying the errors in numerical computation. These are inherent errors, local round-off errors, and local truncation errors.

Inherent Errors

It is that quantity of error which is present in the statement of the problem itself, before finding its solution. It arises due to the simplified assumptions made in the mathematical modeling of a problem. It can also arise when the data is obtained from certain physical measurements of the parameters of the problem.

💡 Why this matters: Inherent errors are unavoidable since they exist before any computation begins — they set a lower bound on the achievable accuracy.

Local Round-off Errors

Every computer has a finite word length and therefore it is possible to store only a fixed number of digits of a given input number. Since computers store information in binary form, storing an exact decimal number in its binary form into the computer memory gives an error. This error is computer dependent.

At the end of computation of a particular problem, the final results in the computer, which is obviously in binary form, should be converted into decimal form — a form understandable to the user — before their print out. Therefore, an additional error is committed at this stage too. This error is called local round-off error.

🔑 Definition — Error: Error = True value – Computed value

🔑 Definition — Absolute error: Absolute error = |Error|

🔑 Definition — Relative error: Relative error = Error / True value

📌 Example: The decimal number (0.7625)₁₀ equals (0.110000110011)₂ in binary. If a computer system has a word length of 12 bits only, then the decimal number 0.7625 is stored in the computer memory in binary form as 0.110000110011. However, it is equivalent to 0.76245. Thus, in storing the number 0.7625, we have committed an error equal to 0.00005, which is the round-off error; inherent with the computer system considered.

Local Truncation Error

It is generally easier to expand a function into a power series using Taylor series expansion and evaluate it by retaining the first few terms. For example, we may approximate the function f(x) = cos x by the series: cos x = 1 − x²/2! + x⁴/4! − ⋯ + (−1)ⁿ x²ⁿ/(2n)! + ⋯

If we use only the first three terms to compute cos x for a given x, we get an approximate answer. Here, the error is due to truncating the series. Suppose, we retain the first n terms, the truncation error (TE) is given by: TE ≤ x²ⁿ⁺² / (2n + 2)!

The TE is independent of the computer used.

If we wish to compute cos x accurate with five significant digits, the question is, how many terms in the expansion are to be included? In this situation: x²ⁿ⁺² / (2n + 2)! < 0.5 × 10⁻⁵ = 5 × 10⁻⁶

Taking logarithm on both sides, we get: (2n + 2) log x − log[(2n + 2)!] < log₁₀ 5 − 6 log₁₀ 10 = 0.699 − 6 = −5.3 or log[(2n + 2)!] − (2n + 2) log x > 5.3

We can observe that the above inequality is satisfied for n = 7. Hence, seven terms in the expansion are required to get the value of cos x, with the prescribed accuracy.

The truncation error is given by: TE ≤ x¹⁶ / 16!

📌 Example: To compute cos x with five significant digits accuracy, one must retain seven terms in the Taylor series expansion because the inequality log[(2n+2)!] − (2n+2) log x > 5.3 is satisfied only when n=7.

💡 Why this matters: Truncation error is analyst-controlled — by choosing how many terms to include, we can achieve any desired accuracy, though at the cost of more computation.

⭐ Key Takeaways

The three fundamental error types in numerical computation are inherent errors (present in the problem statement or data), local round-off errors (due to finite computer word length and binary-decimal conversion), and local truncation errors (due to approximating infinite processes by finite ones). The error is defined as True value minus Computed value, with absolute error being its magnitude and relative error being the error divided by the true value. Round-off errors are computer-dependent, while truncation errors are independent of the computer and can be controlled by retaining more terms. To achieve five significant digits of accuracy when computing cos x using its Taylor series, seven terms are required, as determined by solving the inequality involving the truncation error bound.

🧠 Quick Revision Questions

1. What are the three types of errors discussed in numerical computations, and what is the source of each?

2. Define error, absolute error, and relative error mathematically.

3. Why does the decimal number 0.7625 become 0.76245 in a 12-bit computer system? What type of error does this represent?

4. What is the bound on the truncation error when approximating cos x by the first n terms of its Taylor series?

5. How many terms of the Taylor series for cos x are needed to guarantee five significant digits of accuracy? Show the reasoning.

📘 Lecture 3 — Polynomial

📖 Overview: This lecture introduces the fundamental concepts of polynomial and algebraic equations, along with transcendental equations. It establishes the foundational definitions and key facts necessary for understanding root-finding methods in numerical analysis, explaining why these classification and properties matter for solving equations.

🗂️ Topics Covered

The lecture covers the definition of a polynomial expression of nth degree, the concept of an algebraic equation f(x)=0, and key facts about the roots of such equations—including the number and nature of roots based on degree and sign of the last term. It then defines transcendental equations as those involving logarithmic, trigonometric, or exponential functions, and finally defines the root of an equation as a value 'a' such that f(a)=0.

📝 Lecture Summary

Polynomial

An expression of the form ( f(x) = a_0 x^n + a_1 x^{n-1} + a_2 x^{n-2} + ... + a_{n-1} x + a_n ), where n is a positive integer and ( a_0, a_1, a_2, ..., a_n ) are real constants, is called an nth degree polynomial in x if ( a_0 \neq 0 ). This is the basic building block for algebraic equations in numerical analysis.

🔑 Definition — Polynomial: An expression of the form ( f(x) = a_0 x^n + a_1 x^{n-1} + ... + a_n ) with ( a_0 \neq 0 ), where n is a positive integer and all coefficients are real constants.

Algebraic equation

An equation f(x)=0 is said to be an algebraic equation in x if it is purely a polynomial in x. For example, ( x^5 + x^4 + 3x^2 + x - 6 = 0 ) is a fifth order polynomial and thus an algebraic equation. Other examples include ( x^3 - 6 = 0 ), ( x^6 - 7x = 0 ), ( y^4 - 4y^3 + 3y^2 - y - 2 = 0 ) (polynomial in y), and ( t^4 - 6t^2 - 21 = 0 ) (polynomial in t).

Some facts

1. Every equation of the form f(x)=0 has at least one root; it may be real or complex.
2. Every polynomial of nth degree has n and only n roots.
3. If f(x)=0 is an equation of odd degree, then it has at least one real root whose sign is opposite to that of the last term.
4. If f(x)=0 is an equation of even degree whose last term is negative, then it has at least one positive and at least one negative root.

💡 Why this matters: These facts give a quick way to predict the existence and sign of real roots without solving the equation, which is essential for selecting appropriate numerical methods.

Transcendental equation

An equation is said to be a transcendental equation if it involves logarithmic, trigonometric, and exponential functions, or a combination of all three. For example, ( e^x - 5x - 3 = 0 ), ( e^x - \sin x = 0 ), ( \ln x - \sin x = 0 ), and ( 2 \sec^2 x - \tan x - e^x = 0 ) are all transcendental equations. These are more complex than algebraic equations and often require iterative numerical methods for solution.

🔑 Definition — Transcendental equation: An equation that contains logarithmic, trigonometric, or exponential functions (or any combination thereof).

Root of an equation

For an equation f(x)=0, to find the solution we find such a value which satisfies f(x)=0; these values are known as the roots of the equation. A value 'a' is known as the root of an equation f(x)=0 if and only if f(a)=0.

🔑 Definition — Root of an equation: A value 'a' such that f(a)=0; it satisfies the equation f(x)=0.

📌 Example: For the polynomial ( x^3 - 6 = 0 ), if we find a number 'a' such that ( a^3 - 6 = 0 ), then 'a' is a root of this algebraic equation.

⭐ Key Takeaways

You must distinguish between algebraic (purely polynomial) and transcendental (involving log, trig, or exponential functions) equations. Every nth degree polynomial has exactly n roots (counting multiplicity and complex roots). The sign of the last term and the degree (odd/even) provide critical information about the existence and sign of real roots—odd degree equations always have at least one real root, and even degree equations with a negative last term have both positive and negative real roots. Finally, the root of any equation f(x)=0 is defined precisely by the condition f(a)=0, which is the foundation for all root-finding methods.

🧠 Quick Revision Questions

1. What is the definition of an nth degree polynomial in x?
2. State the four facts about the roots of polynomial equations.
3. What distinguishes a transcendental equation from an algebraic equation?
4. Give three examples of transcendental equations.
5. What is a root of an equation f(x)=0?

📘 Lecture 4 — Properties of an Algebraic Equation

📖 Overview: This lecture introduces fundamental properties of algebraic equations that are essential for root-finding in numerical analysis. It covers the nature of roots, sign rules for determining root signs, and the classification of numerical methods into direct and iterative approaches.

🗂️ Topics Covered

The lecture covers properties of algebraic equations including complex roots occurring in conjugate pairs and factor properties. It explains Descartes' rule of signs for determining positive and negative roots, the intermediate value property for root existence in intervals, and classifies numerical methods into direct methods (like Graefee method) and iterative methods.

📝 Lecture Summary

Properties of an Algebraic equation

An algebraic equation of degree n has specific properties regarding its roots. First, complex roots occur in conjugate pairs — if (a+ib) is a root of f(x)=0, then (a-ib) is also a root. Second, if x=a is a root of the nth-degree polynomial f(x)=0, then (x-a) is a factor of f(x), and dividing f(x) by (x-a) yields a polynomial of degree n-1.

🔑 Definition — Complex conjugate roots: If a polynomial equation has real coefficients and a complex root (a+ib), its complex conjugate (a-ib) is also a root.

Descartes rule of signs

Descartes' rule of signs establishes the relationship between the signs of coefficients and the roots of an equation. The rule states: "The number of positive roots of an algebraic equation f(x)=0 with real coefficients cannot exceed the number of changes in the signs of the coefficients in f(x)=0. Similarly, the number of negative roots cannot exceed the number of changes in the sign of coefficients of f(-x)=0."

🔑 Definition — Descartes' rule of signs: A rule that determines the maximum possible number of positive and negative real roots of a polynomial by counting sign changes in coefficients.

📌 Example: Consider the equation x³ − 3x² + 4x − 5 = 0. In f(x), the coefficient signs are +, −, +, − — there are three sign changes (first + to −, second − to +, third + to −), so at most three positive roots exist. For f(−x) = −x³ − 3x² − 4x − 5, the signs are −, −, −, − — no sign change, so there are no negative roots.

Intermediate value property

The intermediate value property states: If f(x) is a real-valued continuous function in the closed interval a ≤ x ≤ b, and if f(a) and f(b) have opposite signs, then f(x)=0 has at least one root β such that a ≤ β ≤ b. Simply put, if f(a) and f(b) have different signs, the equation must have at least one real root between a and b.

🔑 Definition — Intermediate value property: For a continuous function, if function values at interval endpoints have opposite signs, at least one root exists within that interval.

💡 Why this matters: This property provides the theoretical foundation for bracketing methods like the bisection method, guaranteeing root existence before starting numerical computation.

Numerical methods for solving equations

Numerical methods for solving algebraic or transcendental equations are classified into two groups: Direct methods and Iterative methods.

Direct methods are those that do not require any information about the initial approximation of the root to start the solution. Examples include Graefee root squaring method, Gauss elimination method, and Gauss Jordan method. All these methods do not require any type of initial approximation.

Iterative methods require an initial approximation to start the solution process.

🔑 Definition — Direct methods: Numerical methods that find the solution without needing an initial guess for the root.

🔑 Definition — Iterative methods: Numerical methods that begin with an initial approximation and progressively refine it to approach the actual root.

⭐ Key Takeaways

Complex roots of polynomial equations always occur in conjugate pairs, and if a value is a root, the corresponding linear factor divides the polynomial. Descartes' rule of signs provides a quick way to estimate the maximum number of positive and negative real roots by counting sign changes in f(x) and f(-x). The intermediate value property guarantees at least one root exists in an interval where the function changes sign, forming the basis for bracketing methods. Numerical methods split into direct methods (no initial guess needed, like Graefee and Gauss methods) and iterative methods (require initial approximation). Understanding these properties is essential before applying any root-finding algorithm.

🧠 Quick Revision Questions

1. If a polynomial equation has a complex root (3+2i), what other root must it have?
2. How many sign changes occur in the polynomial x⁴ − 2x³ + 3x² − x + 1 = 0, and what does this indicate about positive roots?
3. For the polynomial f(x) = x³ + 2x² + x + 1, use Descartes' rule to determine the maximum number of negative roots.
4. If f(2) = 4 and f(5) = −3, what does the intermediate value property tell us about roots in the interval [2,5]?
5. Name two direct methods for solving algebraic equations and state what characteristic they share.

📘 Lecture 5 — Bisection Method, Initial Approximations, and Iterative Methods

📖 Overview: This lecture introduces the concept of iterative methods for finding roots of equations and focuses on obtaining initial approximations through graphical and analytical methods. It then provides a detailed explanation and application of the Bisection Method (Bolzano), including procedural steps, examples, and criteria for stopping iterations.

🗂️ Topics Covered

The lecture covers iterative methods, methods for obtaining initial approximations (graphical and analytical), the Bisection method (Bolzano) with its procedure and step-by-step examples, and the criteria for stopping the iterative process to achieve a desired accuracy.

📝 Lecture Summary

Bisection method, Newton raphson method, secant method, jacobi method are all examples of iterative methods. How to get an initial approximation?

The initial approximation for finding a root can be found using two methods: the graphical method and the analytical method.

🔑 Definition — Graphical Method: The equation f(x)=0 can be rewritten as f1(x)=f2(x). The initial approximation is taken as the abscissa (x-coordinate) of the point of intersection of the graphs of y=f1(x) and y=f2(x).

📌 Example: For the equation f(x)=x – sin x – 1 = 0, it is rewritten as y=x-1 and y=sin x. The graphs of these two functions intersect at approximately x=1.9, so the initial approximation is 1.9.

🔑 Definition — Analytical Method: This method is based on the intermediate value property. It involves locating two values, a and b, such that f(a) and f(b) have opposite signs. The root is then known to lie between these two points.

📌 Example: Consider the equation f(x)=3x – 1 + sin x = 0. f(0) = -1 and f(1) = 3 – 1 + sin(1 × 180/π) = 2 + 0.84147 = 1.64299. Since f(0) and f(1) have opposite signs, one root lies between 0 and 1.

💡 Why this matters: The analytical method ensures that a root exists in a given interval by checking the sign change of the function at its endpoints.

Bisection method (Bolzano)

The Bisection method is used to locate the root of an equation f(x)=0 within an interval (x0, x1) where f(x0) and f(x1) have opposite signs, i.e., f(x0)f(x1) < 0. This sign change indicates that the graph crosses the x-axis, guaranteeing at least one root in the interval.

The procedure starts by defining the midpoint x2 = (x0 + x1)/2. If f(x2) = 0, then x2 is the root. Otherwise, the root lies either between x0 and x2 or x1 and x2, depending on which pair produces a sign change. The process of halving the interval is repeated until the root is found to the desired accuracy.

📐 Formula: Midpoint: x_new = (x_left + x_right) / 2 → The next approximation is the average of the two interval endpoints that bracket the root.

📌 Example: Solve x³ – 9x + 1 = 0 for the root between x=2 and x=4.

• f(2) = 8 – 18 + 1 = -9; f(4) = 64 – 36 + 1 = 29. Since f(2)f(4) < 0, root lies between 2 and 4.
• x2 = (2+4)/2 = 3; f(3) = 27 – 27 + 1 = 1. Since f(2)f(3) < 0, root lies between 2 and 3.
• x3 = (2+3)/2 = 2.5; f(2.5) = 15.625 – 22.5 + 1 = -5.875.
• x4 = (2.5+3)/2 = 2.75; f(2.75) = -2.9531.
• x5 = (2.75+3)/2 = 2.875; f(2.875) = -1.1113.
• x6 = (2.875+3)/2 = 2.9375; f(2.9375) = -0.0901.
• The process is continued until the desired accuracy is obtained.

📌 Example: Carry out five iterations for f(x)=2x cos(2x) – (x+1)².

• f(-1) = 0.8322 > 0, f(0) = -1 < 0. Root lies between -1 and 0.
• x2 = (0-1)/2 = -0.5
• f(-0.5) = -0.7903 < 0. Root lies between -1 and -0.5.
• x3 = (-0.5-1)/2 = -0.75
• f(-0.75) = -0.1686 < 0. Root lies between -1 and -0.75.
• x4 = (-0.75-1)/2 = -0.875
• f(-0.875) = 0.296275 > 0. Root lies between -0.875 and -0.75.
• x5 = (-0.75-0.875)/2 = -0.8125
• f(-0.8125) = 0.052970 > 0. Root lies between -0.8125 and -0.75.
• x6 = (-0.75-0.8125)/2 = -0.78125

📌 Example: Carry out the first five iterations for f(x)=x cosx – 2x² + 3x – 1, for 1.2 ≤ x ≤ 1.3.

• f(1.2) = 0.1548 > 0, f(1.3) = -0.1322 < 0. Root lies between 1.2 and 1.3.
• x2 = (1.2+1.3)/2 = 1.25; f(1.25) = 0.0191 > 0. Root lies between 1.25 and 1.3.
• x3 = (1.25+1.3)/2 = 1.275; f(1.275) = -0.0545 < 0. Root lies between 1.25 and 1.275.
• x4 = (1.25+1.275)/2 = 1.2625; f(1.2625) = -0.0172 < 0. Root lies between 1.25 and 1.2625.
• x5 = (1.25+1.2625)/2 = 1.25625; f(1.25625) = 0.00108 > 0. Root lies between 1.25625 and 1.2625.
• x6 = (1.25625+1.2625)/2 = 1.259375

When to stop the process of iteration?

The iteration process stops when a specified condition is met. This condition can be either a fixed number of iterations or a desired accuracy.

🔑 Definition — Accuracy: To achieve an accuracy of 10⁻³, the absolute difference between two consecutive approximations is calculated. If this difference is less than or equal to the required accuracy, the process stops.

📌 Example: If two consecutive roots are 2.135769 and 2.135648, the difference is 2.135769 – 2.135648 = 0.000121. Since 0.000121 is less than 0.001 (10⁻³), the accuracy of 10⁻³ is achieved and the process stops.

⭐ Key Takeaways

The most critical points to remember from this lecture are that iterative methods, such as the Bisection method, require an initial approximation, which can be found graphically or analytically by ensuring function values at interval endpoints have opposite signs. The core procedure of the Bisection method is to repeatedly halve the interval bracketing the root by calculating the midpoint and selecting the subinterval where a sign change occurs. All calculations for trigonometric functions in numerical methods are performed in radians mode with π approximated as 3.14. Finally, the iteration process is terminated either after a predetermined number of steps or when the absolute difference between two successive approximations is less than a specified tolerance or accuracy.

🧠 Quick Revision Questions

1. What are the two methods to obtain an initial approximation for a root?
2. State the fundamental condition that must be satisfied for the analytical method to guarantee a root in an interval.
3. In the Bisection method, if f(x₀) and f(x₁) have opposite signs, what is the formula for the next approximation, x₂?
4. In which mode (degrees or radians) are all calculations in numerical analysis performed?
5. How do you determine when to stop the Bisection method if you are asked to find the root with an accuracy of 10⁻⁴?

📘 Lecture 6 — Regula-Falsi Method (Method of False Position) & Method of Iteration

📖 Overview: This lecture covers two root-finding techniques: the Regula-Falsi (false position) method, which uses a chord between two points with opposite signs to approximate a root, and the Iteration method, which rewrites the equation as (x = \phi(x)) for successive approximations. These methods are essential for solving algebraic and transcendental equations numerically.

🗂️ Topics Covered

Regula-Falsi method derivation and formula; convergence observation about initial interval size; worked examples on cubic, transcendental, and mixed exponential-trigonometric equations; Method of Iteration definition; convergence theorem; examples using cubic and trigonometric functions with high accuracy.

📝 Lecture Summary

Regula-Falsi method (Method of false position)

Here we choose two points (x_n) and (x_{n-1}) such that (f(x_n)) and (f(x_{n-1})) have opposite signs. Intermediate value property suggests that the graph of (y=f(x)) crosses the x-axis between these two points and therefore, a root lies between these two points. Thus to find the real root of (f(x)=0) using Regula-Falsi method, we replace the part of the curve between the points (A(x_n, f(x_n))) and (B(x_{n-1}, f(x_{n-1}))) by a chord in the interval and we take the point of intersection of this chord with x-axis as initial approximation.

The equation of the chord joining points A and B is: [\frac{y - f(x_n)}{x - x_n} = \frac{f(x_{n-1}) - f(x_n)}{x_{n-1} - x_n}]

Setting (y=0) gives the first approximation to the root: [x_{n+1} = x_n - \frac{x_n - x_{n-1}}{f(x_n) - f(x_{n-1})} f(x_n)]

We observe that (f(x_{n-1})) and (f(x_{n+1})) are of opposite signs, thus it is possible to apply the above procedure to determine the line through B and A₁ and so on.

🔑 Definition — Regula-Falsi formula: (x_{n+1} = x_n - \frac{x_n - x_{n-1}}{f(x_n) - f(x_{n-1})} f(x_n)) → finds the x-intercept of the chord joining two points with opposite function signs.

📌 Example: Use Regula-Falsi to compute a real root of (x^3 - 9x + 1 = 0): (i) between 2 and 4, (ii) between 2 and 3.

Let (f(x) = x^3 - 9x + 1).
(f(2) = 8 - 18 + 1 = -9), (f(4) = 64 - 36 + 1 = 29). Opposite signs, root lies between 2 and 4.

First approximation ((x_1=2, x_2=4)): [x_3 = 4 - \frac{4-2}{29-(-9)}(29) = 4 - \frac{2(29)}{38} = 4 - 1.5263 = 2.4736] (f(x_3) = (2.4736)^3 - 9(2.4736) + 1 = 15.13520 - 22.2624 + 1 = -6.12644)

Second approximation: [x_4 = 2.4736 - \frac{2.4736-4}{-6.12644-29}(-6.12644) = 2.4736 + 0.26619 = 2.73989] (f(x_4) = 20.5683 - 24.65901 + 1 = -3.090707)

Third approximation: [x_5 = 2.73989 - \frac{2.73989-4}{-3.090707-29}(-3.090707) = 2.86125] (f(x_5) = 23.42434 - 25.75125 + 1 = -1.326868)

(ii) Between 2 and 3:
(f(2) = -9), (f(3) = 27 - 27 + 1 = 1). Opposite signs, root lies between 2 and 3.

First approximation ((x_1=2, x_2=3)): [x_3 = 3 - \frac{3-2}{1-(-9)}(1) = 3 - \frac{1}{10} = 2.9] (f(x_3) = 24.389 - 26.1 + 1 = -0.711)

Second approximation: [x_4 = 2.9 - \frac{2.9-3}{-0.711-1}(-0.711) = 2.9 + 0.04156 = 2.94156] (f(x_4) = 25.45265 - 26.47404 + 1 = -0.0207)

Third approximation: [x_5 = 2.94156 - \frac{2.94156-3}{-0.0207-1}(-0.0207) = 2.94275] (f(x_5) = 25.48356 - 26.48475 + 1 = -0.0011896)

💡 Why this matters: The value when the interval is (2,3) converges faster and is closer to the root. Important observation: The initial interval ((x_1, x_2)) in which the root lies should be sufficiently small.

Example — Transcendental Equation

Use Regula-Falsi to find a real root of (\ln x - \cos x = 0) accurate to four decimal places after three successive approximations. (Note: all calculations in radians mode.)

(f(x) = \ln x - \cos x)
(f(1) = 0 - 0.540302 = -0.540302 < 0)
(f(2) = 0.69315 - 0.41615 = 1.109)
(f(1)f(2) < 0), root lies between 1 and 2.

First approximation ((x_1=1, x_2=2)): [x_3 = 2 - \frac{2-1}{1.109 - (-0.540302)}(1.109) = 2 - \frac{1.1093}{1.6496} = 1.3275] (f(x_3) = \ln 1.3275 - \cos 1.3275 = 0.2833 - 0.2409 = 0.0424)

Second approximation (using (x_1) and (x_3)): [x_4 = 1.3275 - \frac{(0.3275)(0.0424)}{0.0424 + 0.5403} = 1.3037] (f(x_4) = 1.24816 \times 10^{-3})

Third approximation: [x_5 = 1.3037 - \frac{(0.3037)(0.001248)}{0.001248 + 0.5403} = 1.3030] (f(x_5) = 0.6245 \times 10^{-4})
The required real root is 1.3030.

Example — Mixed Equation

Use method of false position to solve (e^{-x} + 2^{-x} + 2\cos x - 6 = 0) for (1 \le x \le 2).

(f(x) = e^{-x} + 2^{-x} + 2\cos x - 6) with (x_0=1, x_1=2).

(f(1) = 2.7182 + 0.5 + 2(0.5403) - 6 = -1.7011)
(f(2) = 7.3886 + 0.25 + 2(-0.4161) - 6 = 0.8068)

For (n=1): [x_2 = 2 - \frac{2-1}{0.8068 + 1.7011}(0.8068) = 2 - \frac{1}{2.5079}(0.8068) = 1.6783] (f(1.6783) = -0.5457)

For (n=2): [x_3 = 1.6783 - \frac{1.6783-2}{-0.5457-0.8068}(-0.5457) = 1.6783 + 0.12979 = 1.8081] (f(1.8081) = -0.08575)

For (n=3): [x_4 = 1.8081 - \frac{1.8081-1.6783}{-0.08575+0.5457}(-0.08575) = 1.8323] (f(1.8323) = 0.01199)

For (n=4): [x_5 = 1.8323 - \frac{1.8323-1.8081}{0.01199+0.08575}(0.01199) = 1.8323 - 0.00296 = 1.8293] (f(1.8293) = -0.000343)

For (n=5): [x_6 = 1.8293 - \frac{1.8293-1.8323}{-0.000343-0.01199}(-0.000343) = 1.8293]

Example — Trigonometric Equation

Solve (2x\cos 2x - (x-2)^2 = 0) for (2 \le x \le 3), perform three iterations.

(f(x) = 2x\cos 2x - (x-2)^2) with (x_0=2, x_1=3).
(f(2) = 4\cos 4 - 0 = -2.6146)
(f(3) = 6\cos 6 - 1 = 4.7610)

For (n=1): [x_2 = 3 - \frac{3-2}{4.7610 - (-2.6146)}(4.7610) = 3 - \frac{4.7610}{7.3756} = 3 - 0.6455 = 2.3545] (f(2.3545) = 4.709\cos 4.709 - 0.1257 = -0.1416)

For (n=2): [x_3 = 2.3545 - \frac{2.3545-3}{-0.1416-4.7610}(-0.1416) = 2.3731] (f(2.3731) = 4.7462\cos 4.7462 - 0.1392 = -0.1392)

For (n=3): [x_4 = 2.3731 - \frac{2.3731-2.3545}{-0.1392+0.1416}(-0.1392) = 2.3707] (f(2.3707) = 4.7414\cos 4.7412 - 0.1374 = 0.00013)

For (n=4): [x_5 = 2.3707 - \frac{2.3707-2.3731}{0.00013-0.0212}(0.00013) = 2.3707]

Example — Logarithmic Equation

Using Regula-Falsi, find the real root of (x\log_{10}x = 1.2) correct to three decimal places.

(f(x) = x\log_{10}x - 1.2)
(f(2) = 2(0.3010) - 1.2 = -0.5979)
(f(3) = 3(0.4771) - 1.2 = 0.2314)
Opposite signs, root lies between (x_1=2, x_2=3).

First approximation: [x_3 = 3 - \frac{3-2}{0.2314+0.5979}(0.2314) = 3 - \frac{0.2314}{0.8293} = 2.72097] (f(x_3) = 2.72097\log_{10} 2.72097 - 1.2 = -0.01713)

Second approximation: [x_4 = 2.72097 - \frac{2.72097-3}{-0.01713-0.2314}(-0.01713) = 2.7402] (f(x_4) = 2.7402\log_{10} 2.7402 - 1.2 = -3.8905 \times 10^{-4})

Thus, the root correct to three decimal places is 2.740.

💡 Why this matters: If TOL is given, we can find its value by subtracting consecutive roots and writing in exponential notation. Stop when required TOL is obtained.

Method of Iteration

Method of iterations can be applied to find a real root of (f(x)=0) by rewriting it in the form: [x = \phi(x)]

Let (x = x_0) be the initial approximation to the actual root (\alpha). Then: [x_1 = \phi(x_0), \quad x_2 = \phi(x_1), \quad x_3 = \phi(x_2), \quad \ldots, \quad x_n = \phi(x_{n-1})] If the sequence of approximate roots converges to (\alpha), it is taken as the root of (f(x)=0).

🔑 Definition — Iteration method: (x_{n+1} = \phi(x_n)) where (x = \phi(x)) is obtained by rewriting (f(x)=0).

Theorem

If (\alpha) is a root of (f(x)=0) which is equivalent to (x = \phi(x)), and I is any interval containing (x = \alpha) with (|\phi'(x)| < 1) for all (x \in I), then the sequence of approximations (x_1, x_2, x_3, \ldots, x_n) will converge to the root (\alpha) provided that the initial approximation (x_0) is chosen in I.

Example — Trigonometric Iteration

(f(x) = \cos x - 2x + 3 = 0) rewritten as: [x = \frac{1}{2}(\cos x + 3) = \phi(x)] (\phi'(x) = -\frac{1}{2}\sin x)
(f(1) = 1.54030 > 0), (f(2) = -1.41614 < 0), root lies between 1 and 2.
(|\phi'(1)|) and (|\phi'(2)| < 1), so iteration method can be applied.

Let (x_0 = 1.5): (x_1 = \frac{1}{2}(\cos 1.5 + 3) = 1.999825)
(x_2 = \frac{1}{2}(\cos 1.999825 + 3) = 1.999695)
(x_3 = \frac{1}{2}(\cos 1.999695 + 3) = 1.999695)

This is the required root correct up to 5 decimal places.

Example — Cubic Iteration

Find the real root of (x^3 + x^2 - 1 = 0) by iteration.

(f(x) = x^3 + x^2 - 1)
(f(0) = -1 < 0), (f(1) = 1 > 0), root lies between 0 and 1.

Rewrite: (x^2(x+1) = 1 \implies x = \frac{1}{\sqrt{x+1}} = \phi(x))
(\phi'(x) = -\frac{1}{2}(x+1)^{-3/2})
(\phi'(0) = -\frac{1}{2} < 1), (\phi'(1) = -\frac{1}{2\sqrt{8}} < 1), so (|\phi'(x)| < 1) for all values in the interval.

Let (x_0 = 0.65): (x_1 = \frac{1}{\sqrt{1.65}} = 0.7784989)
(x_2 = \frac{1}{\sqrt{1.7784989}} = 0.7498479)
(x_3 = \frac{1}{\sqrt{1.7498479}} = 0.7559617)
(x_4 = \frac{1}{\sqrt{1.7559617}} = 0.7546446)
(x_5 = \frac{1}{\sqrt{1.7546446}} = 0.7549278)
(x_6 = \frac{1}{\sqrt{1.7549278}} = 0.7548668)
(x_7 = \frac{1}{\sqrt{1.7548668}} = 0.7548799)
(x_8 = \frac{1}{\sqrt{1.7548799}} = 0.7548771)
(x_9 = \frac{1}{\sqrt{1.7548771}} = 0.7548777)
(x_{10} = \frac{1}{\sqrt{1.7548777}} = 0.7548776)
(x_{11} = \frac{1}{\sqrt{1.7548776}} = 0.7548776)

Hence root is 0.7548776.

Example — Cosine Iteration

Find a real root of (\cos x = 3x - 1) correct to seven decimal places.

(f(x) = \cos x - 3x + 1)
(f(0) = 1 - 0 + 1 = 1 > 0)
(f(\pi/2) = \cos 1.57 - 4.71 + 1 = 0.0007963 - 4.71 + 1 = -3.7092037 < 0)
Root lies between 0 and (\pi/2).

(\phi(x) = \frac{1}{3}(\cos x + 1)), (\phi'(x) = -\frac{1}{3}\sin x)
Since (\sin x) is bounded between -1 and 1, (|\phi'(x)| < 1), so iteration method can be applied.

Let (x_0 = 0.5): (x_1 = \frac{1}{3}[\cos 0.5 + 1] = 0.6258608)
(x_2 = \frac{1}{3}[\cos 0.6258608 + 1] = 0.6034863)
(x_3 = \frac{1}{3}[\cos 0.6034863 + 1] = 0.6077873)
(x_4 = \frac{1}{3}[\cos 0.6077873 + 1] = 0.6069711)
(x_5 = \frac{1}{3}[\cos 0.6069711 + 1] = 0.6071264)
(x_6 = \frac{1}{3}[\cos 0.6071264 + 1] = 0.6070969)
(x_7 = \frac{1}{3}[\cos 0.6070969 + 1] = 0.6071025)
(x_8 = \frac{1}{3}[\cos 0.6071025 + 1] = 0.6071014)
(x_9 = \frac{1}{3}[\cos 0.6071014 + 1] = 0.6071016)
(x_{10} = \frac{1}{3}[\cos 0.6071016 + 1] = 0.6071016)

Hence root is 0.6071016.

⭐ Key Takeaways

The Regula-Falsi method uses a chord between two points with opposite signs to approximate the root, and its formula is derived from the equation of a line. The initial interval must be sufficiently small for faster convergence. The Method of Iteration rewrites (f(x)=0) as (x=\phi(x)) and generates successive approximations; it converges if (|\phi'(x)|<1) in the interval containing the root. Both methods are applied to algebraic, transcendental, and mixed equations, with iteration commonly requiring more steps but offering clear convergence criteria.

🧠 Quick Revision Questions

1. Derive the Regula-Falsi formula from the equation of the chord joining two points on the curve.
2. Why must the initial interval in Regula-Falsi be sufficiently small?
3. State the convergence condition for the Method of Iteration and explain its significance.
4. Rewrite (x^3 - 4x + 1 = 0) in the form (x = \phi(x)) suitable for iteration, and verify the convergence condition near the root in (0,1).
5. How do you determine when to stop iterating in both methods when a tolerance (TOL) is given?

📘 Lecture 7 — Newton-Raphson Method

📖 Overview: This lecture introduces the Newton-Raphson method, a powerful and widely used technique for finding roots of nonlinear equations. It derives the iterative formula from Taylor series, provides a geometric interpretation using tangents, and discusses convergence properties, drawbacks, and practical applications through solved examples.

🗂️ Topics Covered

The lecture covers the derivation of the Newton-Raphson formula from Taylor series expansion, its geometrical interpretation using tangents to the curve, important notes on its application and limitations, step-by-step solved examples for finding roots of polynomial and transcendental equations, drawbacks including divergence at inflection points and division by zero, and a detailed analysis of its quadratic convergence rate.

📝 Lecture Summary

Newton-Raphson Method

This method is used for finding a root of f(x)=0. The formula is derived by using the first two terms in Taylor's series expansion. Setting f(x_{n+1}) = 0 and simplifying gives the Newton-Raphson iterative formula.

🔑 Definition — Newton-Raphson Formula: x_{n+1} = x_n - f(x_n) / f'(x_n) for n=0,1,2,… This gives successive approximations to the root.

Geometrical Interpretation

If the curve f(x)=0 meets the x-axis at x=α (the root), and x_0 is a point near α, then the tangent at P_0[x_0, f(x_0)] is: y - f(x_0) = f'(x_0)(x - x_0). This tangent cuts the x-axis at x_1 = x_0 - f(x_0)/f'(x_0), which is the first approximation. Repeating this process gives better approximations rapidly.

Important Notes

1. When f'(x) is very large (slope is large), h will be small, and the root can be found quickly.
2. If the initial approximation x_0 is close to the root, convergence is very fast.
3. The method fails if f'(x) = 0 in the neighborhood of the root. In such cases, the Regula Falsi method should be used.
4. If the initial approximation is not given, choose two points a and b such that f(a) and f(b) have opposite signs. If |f(a)| < |f(b)|, use a as the initial approximation.
5. Newton-Raphson is also called the method of tangents.

Example 1: Solve x^3 - x - 1 = 0

Find a real root correct to four decimal places.

Given f(x) = x^3 - x - 1

• f(1) = -1 < 0, f(2) = 5 > 0, so root lies between 1 and 2.
• f'(x) = 3x^2 - 1
• f''(x) = 6x
• f(2) and f''(2) have the same sign (+), so x_0 = 2

📌 Iteration 1: x_1 = 2 - 5/11 = 1.54545 📌 Iteration 2: f(1.54545) = 1.14576, f'(1.54545) = 6.16487, x_2 = 1.54545 - 1.14576/6.16525 = 1.35961 📌 Iteration 3: f(1.35961) = 0.15369, f'(1.35961) = 4.54562, x_3 = 1.35961 - 0.15369/4.54562 = 1.32579 📌 Iteration 4: f(1.32579) = 4.60959×10^-3, f'(1.32579) = 4.27316, x_4 = 1.32579 - 0.00460959/4.27316 = 1.32471 📌 Iteration 5: f(1.32471) = -3.39345×10^-5, f'(1.32471) = 4.26457, x_5 = 1.32471 + 0.0000339345/4.26457 = 1.324718 Since f(x_5) = 1.823×10^-7, the required root is 1.3247.

Comparison with Bracketing Methods

Methods like Bisection and False Position require bracketing the root and are always convergent. The Newton-Raphson method is an open method — it needs only one initial guess and is not guaranteed to converge, but when it does, it converges very quickly.

Drawbacks of Newton-Raphson Method

1. Divergence at inflection points: If a guess is near an inflection point (where f''(x) = 0), the iterations may diverge.
2. Division by zero or near-zero: If f'(x_i) is very small, x_{i+1} can become very large.
3. Root jumping: For oscillating functions with multiple roots, an initial guess close to one root may converge to a different root.
4. Oscillations near local max/min: Iterations may oscillate about a local maximum or minimum without converging to a root, eventually causing division by near-zero.

Convergence of Newton-Raphson Method

Comparing the N-R formula with the general iteration formula x_{n+1} = φ(x_n), we have φ(x) = x - f(x)/f'(x). The method converges if |φ'(x)| < 1, which requires f(x) f''(x) < [f'(x)]^2 in the interval of interest.

Definition: Let x_n = α + ε_n and x_{n+1} = α + ε_{n+1}, where α is the root. If ε_{n+1} = K ε_n^p, then p is the rate of convergence.

The N-R method has quadratic convergence (p=2). Using Taylor expansion and the fact that f(α)=0, we derive: ε_{n+1} = [ε_n^2 f''(α)] / [2 f'(α)]

So ε_{n+1} = K ε_n^2 where K = f''(α) / [2 f'(α)]. This shows second-order (quadratic) convergence.

Example: Square Root of a Positive Number N

To find √N, solve f(x) = x^2 - N = 0. With f'(x) = 2x, the N-R iteration is: x_{n+1} = x_n - (x_n^2 - N)/(2x_n) = (1/2)(x_n + N/x_n)

📌 Example: Find √12

• Initial guess: x_0 = (3 + 4)/2 = 3.5
• x_1 = (1/2)(3.5 + 12/3.5) = 3.4643
• x_2 = (1/2)(3.4643 + 12/3.4643) = 3.4641
• x_3 = (1/2)(3.4641 + 12/3.4641) = 3.4641 Hence, √12 = 3.4641

Example: Solve f(x) = x - 0.8 - 0.2 sin x in [0, π/2]

Given f(x) = x - 0.8 - 0.2 sin x

• f(0) = -0.8, f(1.57) ≈ 0.570, so root exists.
• f'(x) = 1 - 0.2 cos x
• Since |f(0)| is less, take x_0 = 1.57? No — the note says if |f(a)| < |f(b)|, take a. Here we start with x_0 = 0 based on |f(0)| = 0.8 vs |f(1.57)|=0.57; 0.57 < 0.8, so x_0 = 1.57? The text uses x_0 = 0 differently. Let’s follow the worked solution.

📌 Iteration 1: x_1 = 0 - (-0.8)/0.8 = 1 📌 Iteration 2: f(1)=0.0317, f'(1)=0.8919, x_2 = 1 - 0.0317/0.8919 = 0.9645 📌 Iteration 3: f(0.9645)=0.0002, f'(0.9645)=0.88604, x_3 = 0.9645 - 0.0002/0.88604 = 0.9643

Example: Solve f(x) = 4x cos x - (x-2)^2 in [0, 8]

Given f(x) = 4x cos x - (x-2)^2

• f(0) = -4, f(8) ≈ -66.645, so both negative — no sign change.
• f'(x) = 4 cos 2x - 8x sin 2x - 2(x-2)
• Since |f(8)| is larger, x_0 = 8

📌 Iteration 1: x_1 = 8 - (-66.645)/2.5952 = 33.6801 📌 Iteration 2: x_2 = 33.6801 + 1028.25/200.79 = 38.8011 📌 Iteration 3: x_3 = 38.8011 + 1446.14/(-326.205) = 43.2343 (Note: This example shows divergence — the iterates are moving away from the [0,8] interval.)

Example: Solve ln(x-1) + cos(x-1) = 0 in [1.2, 2]

Given f(x) = ln(x-1) + cos(x-1)

• f(1.2) = -0.6293, f(2) = 0.5403
• f'(x) = 1/(x-1) - sin(x-1)
• f'(1.2) = 5 - 0.1986 = 4.8014, so x_0 = 1.2

📌 Iteration 1: x_1 = 1.2 - (-0.6293)/4.8014 = 1.3311 📌 Iteration 2: f(1.3311)=-0.1596, f'(1.3311)=2.6951, x_2 = 1.3311 - (-0.1596)/2.6951 = 1.3903 📌 Iteration 3: f(1.3903)=-0.016, f'(1.3903)=2.1816, x_3 = 1.3903 - (-0.016)/2.1816 = 1.3976

💡 Why this matters: The Newton-Raphson method is essential for quickly finding roots of complex equations where analytical solutions are impossible, such as in engineering, physics, and economics. Its quadratic convergence makes it highly efficient when a good initial guess is provided.

⭐ Key Takeaways

The Newton-Raphson method is a powerful open iterative method for finding roots of f(x)=0, using the formula x_{n+1} = x_n - f(x_n)/f'(x_n), derived from Taylor series expansion. It has quadratic convergence, meaning the error reduces as the square of the previous error, making it very fast when the initial guess is close to the root. However, the method has critical limitations: it fails when f'(x)=0 near the root, may diverge at inflection points, and can exhibit root jumping or oscillations near local maxima/minima. For optimal results, the initial approximation should be chosen carefully — if bracketing is possible, choose the endpoint where f(x) and f''(x) have the same sign. The method is particularly useful for extracting square roots and solving transcendental equations.

🧠 Quick Revision Questions

1. Derive the Newton-Raphson formula using the first two terms of Taylor's series expansion. What assumption is made about f(x_{n+1})?

2. Explain the geometrical interpretation of the Newton-Raphson method. How does the tangent line at each iteration lead to the improved estimate?

3. What are the four main drawbacks of the Newton-Raphson method? Under what conditions does each problem occur?

4. Show that Newton-Raphson method has quadratic convergence. Starting from x_n = α + ε_n, derive the relationship ε_{n+1} = K ε_n^2 and find the expression for K.

5. Using Newton's method, derive the iterative scheme for finding the square root of a positive number N. Then find √12 correct to four decimal places using this scheme.

📘 Lecture 8 — Secant Method

📖 Overview: This lecture covers the Secant Method, a modification of the Newton-Raphson method for finding roots of equations. It explains the derivation, geometrical interpretation, convergence properties, and provides worked examples. The lecture also introduces Muller’s Method, which uses quadratic interpolation to approximate roots, including complex ones.

🗂️ Topics Covered

The lecture begins with the Secant Method, its derivation by replacing the derivative in Newton-Raphson with a difference ratio, and its geometrical interpretation as drawing secants. It then discusses the method's convergence rate (super-linear, order ~1.618) and provides a detailed step-by-step example for finding the root of a cubic equation. A second example for a transcendental equation is also worked out. Finally, the lecture introduces Muller’s Method, which approximates the function f(x)=0 by a second-degree polynomial (quadratic) passing through three points, allowing it to find both real and complex roots, and derives the formula for the next approximation.

📝 Lecture Summary

Secant Method

The Secant Method is a modified form of the Newton-Raphson method. In the Newton-Raphson method, the derivative f'(xₙ) is replaced by a difference ratio:

f'(xₙ) = [f(xₙ) – f(xₙ₋₁)] / (xₙ – xₙ₋₁)

Where xₙ and xₙ₋₁ are two approximations of the root. Substituting this into the Newton-Raphson formula gives the secant formula:

xₙ₊₁ = xₙ – [f(xₙ)(xₙ – xₙ₋₁)] / [f(xₙ) – f(xₙ₋₁)] = [xₙ₋₁ f(xₙ) – xₙ f(xₙ₋₁)] / [f(xₙ) – f(xₙ₋₁)]

provided f(xₙ) ≠ f(xₙ₋₁). This method requires two starting values, x₀ and x₁, and the function values at both these points are calculated. The next point on the curve is obtained by using the derived formula. This procedure is continued until the root is found with the required accuracy.

📌 Example: Do three iterations of secant method to find the root of f(x) = x³ – 3x + 1 = 0, taking x₀ = 1, x₁ = 0.5.

• n = 1: f(x₀) = f(1) = 1³ – 3(1) + 1 = -1; f(x₁) = f(0.5) = 0.5³ – 3(0.5) + 1 = -0.375. x₂ = [x₀ f(x₁) – x₁ f(x₀)] / [f(x₁) – f(x₀)] = [(1)(-0.375) – (0.5)(-1)] / [-0.375 – (-1)] = 0.125 / 0.625 = 0.2
• n = 2: f(x₂) = f(0.2) = 0.2³ – 3(0.2) + 1 = 0.408. x₃ = [x₁ f(x₂) – x₂ f(x₁)] / [f(x₂) – f(x₁)] = [(0.5)(0.408) – 0.2(-0.375)] / [0.408 – (-0.375)] = 0.279 / 0.783 = 0.3563
• n = 3: f(x₃) = f(0.3563) = 0.3563³ – 3(0.3563) + 1 = -0.02367. x₄ = [x₂ f(x₃) – x₃ f(x₂)] / [f(x₃) – f(x₂)] = [(0.2)(-0.02367) – 0.3563(0.408)] / [-0.02367 – 0.408] = -0.1498 / -0.43167 = 0.3471 x₅ = 0.3473, f(x₅) = -0.0000096, and x₅ - x₄ = 0.0004. Though X5 is not the true root, it is a good approximation and convergence is faster than bisection.

💡 Why this matters: The Secant method converges faster than linear (Bisection) and slower than Newton’s quadratic rate, but it avoids the need to compute the derivative analytically.

Geometrical Interpretation

Geometrically, the secant method corresponds to drawing secants rather than tangents to obtain various approximations to root α. To obtain x₂, we find the intersection between the secant through the points (x₀, f(x₀)) and (x₁, f(x₁)) and the x-axis.

🔑 Definition — Secant Method Geometry: The next approximation x₂ is obtained as the point of intersection of y = 0 and the chord passing through the points (x₀, f(x₀)) and (x₁, f(x₁)).

The equation of the chord is: y – f(x₀) = [f(x₁) – f(x₀)] / (x₁ – x₀) * (x – x₀)

Putting y = 0 gives: x = x₂ = [x₀ f(x₁) – x₁ f(x₀)] / [f(x₁) – f(x₀)]

The method converges rather quickly, though there is a possibility of divergence if the two roots lie on the same side of the curve. The order of convergence of the secant method is (1 + √5)/2 ≈ 1.618. This shows that this method has an order of convergence slightly inferior to that of Newton-Raphson method. In this method, f(x) is not required to change signs between the estimates.

💡 Why this matters: The secant method does not require the function to change sign between successive approximations, which is a requirement in the Regula-Falsi method.

Convergence of Secant Method

The sequence is generated by the rule: xₙ₊₁ = [xₙ₋₁ f(xₙ) – xₙ f(xₙ₋₁)] / [f(xₙ) – f(xₙ₋₁)]

Starting with x₀ and x₁, the sequence {x₀, x₁, ...} will converge to ‘a’, where f(a) = 0.

The Secant method converges faster than linear and slower than Newton’s quadratic convergence.

Muller’s Method

In Muller’s method, f(x) = 0 is approximated by a second-degree polynomial; that is, by a quadratic equation that fits through three points in the vicinity of a root. The roots of this quadratic equation are then approximated to the roots of the equation f(x) = 0. This method is iterative in nature and does not require the evaluation of derivatives as in Newton-Raphson method. This method can also be used to determine both real and complex roots of f(x) = 0.

Suppose xᵢ₋₂, xᵢ₋₁, xᵢ be any three distinct approximations to a root of f(x) = 0, with function values f(xᵢ₋₂) = fᵢ₋₂, f(xᵢ₋₁) = fᵢ₋₁, f(xᵢ) = fᵢ. A general polynomial of second degree is given by f(x) = ax² + bx + c. The coefficients a, b, c are determined by the condition that the polynomial passes through the three points.

We define: hᵢ = xᵢ – xᵢ₋₁ hᵢ₋₁ = xᵢ₋₁ – xᵢ₋₂ λᵢ = hᵢ / hᵢ₋₁ δᵢ = 1 + λᵢ

Using these substitutions, a simplified equation is obtained. To compute λ, set f = 0, which leads to a quadratic equation in λ. To avoid loss of accuracy, a specific formula is used: λ = -2 fᵢ δᵢ / { gᵢ ± [gᵢ² – 4 fᵢ δᵢ λᵢ (fᵢ₋₂ λᵢ – fᵢ₋₁ δᵢ + fᵢ) ]¹/² }

Here, the positive sign must be chosen so that the denominator becomes largest in magnitude. A better approximation to the root is then obtained by using: xᵢ₊₁ = xᵢ + hᵢ λ

⭐ Key Takeaways

The Secant Method is a root-finding technique derived by replacing the derivative in Newton-Raphson with a finite difference ratio, requiring two initial guesses but no derivative evaluation. Its order of convergence is approximately 1.618, which is faster than linear methods like Bisection but slower than Newton's quadratic rate; it does not require the function to change sign between successive approximations. A full step-by-step example for the cubic equation x³ – 3x + 1 = 0 demonstrates the iterative formula xₙ₊₁ = [xₙ₋₁ f(xₙ) – xₙ f(xₙ₋₁)] / [f(xₙ) – f(xₙ₋₁)]. Muller’s Method is an alternative that approximates the function by a quadratic through three points, allowing it to find complex roots without derivative evaluation. The key difference from Regula-Falsi is that the secant method does not check whether the root lies between successive iterates, which can lead to divergence if points are on the same side of the curve.

🧠 Quick Revision Questions

1. Derive the formula for the Secant Method starting from the Newton-Raphson formula.
2. What is the order of convergence of the Secant Method, and how does it compare to the Bisection and Newton-Raphson methods?
3. Perform two iterations of the Secant Method to approximate a root of f(x) = x³ – 2x – 5 = 0, using initial guesses x₀ = 2 and x₁ = 3.
4. What is the key difference between the Secant Method and the Regula-Falsi method in terms of the condition on the root's location between successive approximations?
5. Explain the general principle of Muller’s Method and state one key advantage it has over the Secant Method.

📘 Lecture 9 — Muller’s Method

📖 Overview: This lecture covers two advanced methods for finding roots of equations. Muller's method generalizes the secant method using three starting points and quadratic interpolation to converge faster than the secant method. Graeffe's Root Squaring Method is introduced as a powerful technique for finding all roots of a polynomial simultaneously by repeatedly squaring the polynomial and taking appropriate roots.

🗂️ Topics Covered

Muller's method is presented as an iterative method requiring three starting points, constructing a parabola through them, and using the quadratic formula for the next approximation. A complete two-iteration example is worked through for solving x³ − 3x + 1 = 0. Graeffe's Root Squaring Method is introduced for finding all roots of polynomial equations by repeatedly squaring the polynomial, transforming roots to their powers, and then extracting roots. Multiple examples demonstrate the method including cubic and quartic equations. The lecture concludes with a derivation of the Newton-Raphson extended formula (Chebyshev's formula of third order).

📝 Lecture Summary

Muller’s Method

Muller's method is a generalization of the secant method. It does not require the derivative of the function. It is an iterative method that requires three starting points x₀, x₁, and x₂. A parabola is constructed that passes through the three points; then the quadratic formula is used to find a root of the quadratic for the next approximation. Near a simple root, Muller's method converges faster than the secant method and almost as fast as Newton's method. The method can find real or complex zeros.

💡 Why this matters: Unlike Newton's method, Muller's method does not require derivatives, making it useful when derivatives are difficult or impossible to compute.

Example: Do two iterations of Muller's method to solve x³ − 3x + 1 = 0 starting with x₀ = 0.5, x₁ = 0, x₂ = 1

Solution:

Iteration 1:

• f(x₀) = f(0.5) = (0.5)³ − 3(0.5) + 1 = −0.375
• f(x₁) = f(1) = (1)³ − 3(1) + 1 = −1
• f(x₂) = f(0) = 0 − 3(0) + 1 = +1
• c = f₀ = −0.375
• h₁ = x₁ − x₀ = 0.5 − 0.5 = 0
• h₂ = x₀ − x₂ = 0.5 − 0 = 0.5

Wait — correction from lecture:

• h₁ = x₁ − x₀ = 0.5 (Note: x₁ = 1, x₀ = 0.5, so h₁ = 0.5)
• h₂ = x₀ − x₂ = 0.5 − 0 = 0.5

a = [h₂f₁ − (h₁ + h₂)f₀ + h₁f₂] / [h₁h₂(h₁ + h₂)] = [(0.5)(−1) − (−0.375)(0.5 + 0.5) + (0.5)(1)] / [(0.5)(0.5)(1)] = [−0.5 − (−0.375)(1) + 0.5] / [0.25] = [−0.5 + 0.375 + 0.5] / 0.25 = 0.375 / 0.25 = 1.5

b = (f₁ − f₀ − ah₁²) / h₁ = (−1 − (−0.375) − (1.5)(0.5)²) / 0.5 = (−1 + 0.375 − 0.375) / 0.5 = (−1) / 0.5 = −2

x₃ = x₀ − [2c / (b − √(b² − 4ac))] = 0.5 − [2(−0.375) / (−2 − √(4 − 4(1.5)(−0.375)))] = 0.5 − [−0.75 / (−2 − √(4 + 2.25))] = 0.5 − [−0.75 / (−2 − √6.25)] = 0.5 − [−0.75 / (−2 − 2.5)] = 0.5 − [−0.75 / (−4.5)] = 0.5 − 0.16667 = 0.33333

Iteration 2: Take x₂ = 0, x₀ = 0.33333, x₁ = 0.5

• h₁ = x₁ − x₀ = 0.5 − 0.33333 = 0.16667
• h₂ = x₀ − x₂ = 0.33333 − 0 = 0.33333
• c = f₀ = f(0.33333) = (0.33333)³ − 3(0.33333) + 1 = 0.037046
• f₁ = (0.5)³ − 3(0.5) + 1 = −0.375
• f₂ = 0³ − 3(0) + 1 = 1

a = [(0.33333)(−0.375) − (0.16667 + 0.33333)(0.037046) + (0.16667)(1)] / [(0.16667)(0.33333)(0.5)] = [−0.125 + 0.018523 + 0.16667] / 0.027778 = 0.023148 / 0.027778 = 0.8333

b = (−0.375 − 0.037046 − (0.8333)(0.16667)²) / 0.16667 = −2.6

x₃ = 0.33333 − [2(0.037046) / (−2.6 − √(6.76 − 4(0.8333)(0.037046)))] = 0.33333 − [0.074092 / (−2.6 − √(6.76 − 0.1235))] = 0.33333 − [0.074092 / (−2.6 − √6.6365)] = 0.33333 − [0.074092 / (−2.6 − 2.5761)] = 0.33333 − [0.074092 / (−5.1761)] = 0.33333 + 0.01431 = 0.34764

For third iteration, take x₂ = 0.333333, x₀ = 0.3475, x₁ = 0.5

Graeffe’s Root Squaring Method

Graeffe's Root Squaring Method is particularly attractive for finding all the roots of a polynomial equation. Consider a polynomial of third degree: f(x) = a₀ + a₁x + a₂x² + a₃x³

f(x)f(−x) = a₃²x⁶ − (a₂² − 2a₁a₃)x⁴ + (a₁² − 2a₀a₂)x² − a₀²

With substitution t = x²: f(x)f(−x) = a₃²t³ − (a₂² − 2a₁a₃)t² + (a₁² − 2a₀a₂)t − a₀²

The roots of this equation are squares or 2ⁱ (i = 1) powers of the original roots. After each squaring, the coefficients become large and overflow is possible as i increases. After squaring 'i' times, we estimate the roots by evaluating the 2ⁱ-th root of aᵢ₋₁/aᵢ.

🔑 Definition — Graeffe's Root Squaring Method: A method for finding all roots of a polynomial by repeatedly squaring the polynomial to separate roots by magnitude, then extracting appropriate roots to recover the original roots.

The method fails when the roots of the given polynomial are repeated. It does not require prior information about approximate values of the roots but is applicable to polynomials only.

Procedure: For polynomial f(x) = xⁿ + a₁xⁿ⁻¹ + a₂xⁿ⁻² + ... + aₙ₋₁x + aₙ, separate even and odd powers, square both sides, put x² = y to get new equation yⁿ + b₁yⁿ⁻¹ + ... + bₙ = 0 where b₁ = −a₁² + 2a₂, b₂ = a₂² − 2a₁a₃ + 2a₄, ..., bₙ = (−1)ⁿaₙ².

If p₁, p₂, ..., pₙ are roots of original equation, then roots of the squared equation are p₁², p₂², ..., pₙ².

Example: Using Graeffe root squaring method, find all the roots of x³ − 6x² + 11x − 6 = 0

Solution: For i = 1: x³ − (36 − 22)x² + (121 − 72)x − 36 = x³ − 14x² + 49x − 36 For i = 2: x³ − (196 − 98)x² + (2401 − 1008)x − 1296 = x³ − 98x² + 1393x − 1296 For i = 3: x³ − (9604 − 2786)x² + (1940449 − 254016)x − 1679616 = x³ − 6818x² + 16864333x − 1679616

Root estimates from polynomial (i=1): ⁴√(36/49) = 0.85714, ⁴√(49/14) = 1.8708, ⁴√(14/1) = 3.7417

From polynomial (i=2): ⁸√(1296/1393) = 0.9821, ⁸√(1393/98) = 1.9417, ⁸√(98/1) = 3.1464

From polynomial (i=3): ¹⁶√(1679616/16864333) = 0.99949, ¹⁶√(16864333/6818) = 1.99143, ¹⁶√(6818/1) = 3.0144

Exact roots: 1, 2, 3

Example: Apply Graeffe's root squaring method to solve x³ − 8x² + 17x − 10 = 0

Three sign changes are observed, so f(x) may have three positive roots.

After repeated squaring and substitution, the final equation is u³ − 378882u² + 264082u − 10⁸ = 0

p₁ = (378882)^(1/8) = 4.9809593 ≈ 5 p₂ = (264082/378882)^(1/8) = 0.9558821 ≈ 1 p₃ = (10⁸/264082)^(1/8) = 2.1003064 ≈ 2

Check: f(5) = f(2) = f(1) = 0, confirming all are roots.

Example: Find all roots of x⁴ − 3x + 1 = 0

Two sign changes → two positive real roots. No sign change in f(−x) → two complex roots.

After repeated squaring and substitution, final equation: u⁴ − 4u³ + 645u² − 5917u + 1 = 0

p₁ = 4^(1/8) = 1.189207 (real) p₄ = (1/5971)^(1/8) = 0.3376659 (real)

For complex roots ρe^(±iφ) = ξ ± iη: ρ = (5917/654)^(1/4) = 1.5780749 ξ₂ = −½(p₁ + p₄) = −0.7634365 η₂ = √(ρ² − ξ²) = √(1.9074851) = 1.3811173

The four roots: 1.1892071, 0.3376659, −0.7634365 ± 1.3811173i

Revision Example: Newton-Raphson Extended Formula

Obtain the Newton-Raphson extended formula for finding the root of f(x) = 0:

x₁ = x₀ − f(x₀)/f'(x₀) − ½[f(x₀)]²f''(x₀)/[f'(x₀)]³

Solution: Expanding f(x) by Taylor's series and retaining up to second-order term: 0 = f(x) = f(x₀) + (x − x₀)f'(x₀) + (x − x₀)²f''(x₀)/2

f(x₁) = f(x₀) + (x₁ − x₀)f'(x₀) + (x₁ − x₀)²f''(x₀)/2 = 0

This gives the Newton-Raphson extended formula: x₁ = x₀ − f(x₀)/f'(x₀) − ½[f(x₀)]²f''(x₀)/[f'(x₀)]³

🔑 Definition — Chebyshev's formula of third order: An extended Newton-Raphson formula that includes second derivative information, providing cubic convergence.

⭐ Key Takeaways

Muller's method uses three points to construct a parabola and requires no derivatives, converging faster than secant method but almost as fast as Newton's method near simple roots. Graeffe's Root Squaring Method finds all roots of a polynomial simultaneously by repeatedly squaring the polynomial, separating roots by magnitude, then extracting the 2ⁱ-th root of coefficient ratios. The method handles real and complex roots but fails for repeated roots. The Newton-Raphson extended formula (Chebyshev's formula) improves convergence by incorporating second derivative information, giving cubic instead of quadratic convergence. Graeffe's method is especially useful because it requires no initial approximations, though it only works on polynomials and coefficients may overflow.

🧠 Quick Revision Questions

1. How many starting points does Muller's method require, and what geometric object is constructed through them?
2. In Graeffe's Root Squaring Method, what relationship exists between the roots of the transformed equation and the original equation?
3. What is the convergence rate of Muller's method near a simple root, and how does it compare to Newton's method and the secant method?
4. For the equation x³ − 6x² + 11x − 6 = 0, the Graeffe method with i=3 gives root approximations of 0.99949, 1.99143, and 3.0144. What are the exact roots, and what happens to the approximations as i increases?
5. Write down the Newton-Raphson extended formula (Chebyshev's formula of third order) and explain what makes it different from the standard Newton-Raphson formula.

📘 Lecture 10 — Solution of Linear System of Equations and Matrix Inversion

📖 Overview: This lecture introduces numerical methods for solving systems of linear equations. It covers direct methods (Gaussian elimination with partial pivoting) for obtaining exact solutions after a finite number of steps, explaining the process of reducing a system to upper triangular form and using back substitution. The importance of pivoting to avoid zero divisors and ensure numerical stability is emphasized through detailed examples.

🗂️ Topics Covered

The lecture begins by classifying numerical methods into direct and iterative categories. It then focuses on the Gaussian elimination method, explaining the two stages: reduction to upper triangular form using elementary row transformations, and back substitution. The concept of pivoting, including partial and full pivoting, is introduced to handle cases where pivot elements become zero. Several fully worked examples demonstrate the application of Gaussian elimination with and without partial pivoting to solve 3×3 and 4×4 systems.

📝 Lecture Summary

Solution of Linear System of Equations and Matrix Inversion

A system of n linear equations in n unknowns has a unique solution if the determinant of the coefficient matrix A is non-zero (|A| ≠ 0). If |A| = 0, either no solution or infinite solutions exist. Numerical methods are classified into direct methods (where the solution is obtained after performing all steps) and iterative methods (where an initial approximation is improved step-by-step). Direct methods include elimination methods (Gaussian elimination, Gauss-Jordan) and decomposition methods (Crout’s/Cholesky’s reduction). Iterative methods include Jacobi, Gauss-Seidel, and relaxation methods.

💡 Why this matters: The choice between direct and iterative methods depends on the size and structure of the system. Direct methods are preferred for smaller, dense systems, while iterative methods are more efficient for large, sparse systems.

Gaussian Elimination Method

This method has two stages. Stage I: The given system is reduced to an equivalent upper triangular form using elementary row transformations. Stage II: The upper triangular system is solved using back substitution.

🔑 Definition — Upper triangular form: A matrix where all entries below the main diagonal are zero (cᵢⱼ = 0 for i > j). 📐 Process: Starting with the system Ax = b, the first equation is used to eliminate x₁ from all subsequent equations. This process is repeated for x₂, x₃, etc., until the system takes the form: c₁₁x₁ + c₁₂x₂ + … + c₁ₙxₙ = d₁ c₂₂x₂ + … + c₂ₙxₙ = d₂ ... cₙₙxₙ = dₙ

Then, backward substitution proceeds: xₙ = dₙ/cₙₙ, then substitute to find xₙ₋₁, and so on. 📌 Example: Solve 2x + 3y - z = 5, 4x + 4y - 3z = 3, -2x + 3y - z = 1. Stage I: Divide first eqn by 2, eliminate x from eqn 2 and 3 → x + 1.5y - 0.5z = 2.5, -2y - z = -7, 6y - 2z = 6. Divide 2nd eqn by -2, eliminate y from eqn 3 → x + 1.5y - 0.5z = 2.5, y + 0.5z = 3.5, -5z = -15. Stage II: From eqn 3, z = 3. From eqn 2, y = 3.5 - 0.5(3) = 2. From eqn 1, x = 2.5 - 1.5(2) + 0.5(3) = 1. Solution: x = 1, y = 2, z = 3.

Partial and Full Pivoting

The Gaussian elimination method fails if any pivot element becomes zero. Pivoting is the process of reordering equations (and occasionally variables) to avoid zero pivots.

🔑 Definition — Partial pivoting: If the pivot in the i-th column is zero, search all rows below (j > i) for the numerically largest element in that column, then interchange the i-th row with the j-th row. 🔑 Definition — Full (complete) pivoting: Both rows and columns are interchanged so that the largest element in the entire remaining submatrix becomes the pivot. This is more complicated and less commonly used for hand computation.

📌 Example (Partial pivoting necessity): Solve the system 10⁻⁵x₁ + x₂ = 1, x₁ + x₂ = 2. Without pivoting (using first eqn to eliminate x₁): The multiplier is 10⁵, leading to severe round-off error and an incorrect solution x₁ = 0, x₂ = 1. With partial pivoting: Swap rows so the largest element (1) in column 1 is the pivot → x₁ + x₂ = 2, 10⁻⁵x₁ + x₂ = 1. Now elimination gives the correct solution x₁ = 1, x₂ = 1.

📌 Example: Solve x + y + z = 7, 3x + 3y + 4z = 24, 2x + y + 3z = 16 using Gaussian elimination with partial pivoting. Matrix form: [1 1 1 | 7], [3 3 4 | 24], [2 1 3 | 16]. Pivot check: a₁₁ = 1 (non-zero), but column 1 has a larger element 3 in row 2. Swap R₁ and R₂ → [3 3 4 | 24], [1 1 1 | 7], [2 1 3 | 16]. Stage I: Divide R₁ by 3 → [1 1 4/3 | 8], [1 1 1 | 7], [2 1 3 | 16]. Eliminate x from R₂ (R₂ - R₁) and R₃ (R₃ - 2R₁) → [1 1 4/3 | 8], [0 0 -1/3 | -1], [0 -1 1/3 | 0]. Swap R₂ and R₃ to get non-zero pivot → [1 1 4/3 | 8], [0 -1 1/3 | 0], [0 0 -1/3 | -1]. Stage II: From eqn 3, z = 3. From eqn 2, -y + 1 = 0 → y = 1. From eqn 1, x + 1 + 4 = 8 → x = 3.

📌 Example: Solve 0x₁ + 4x₂ + 2x₃ + 8x₄ = 24, 4x₁ + 10x₂ + 5x₃ + 4x₄ = 32, 4x₁ + 5x₂ + 6.5x₃ + 2x₄ = 26, 9x₁ + 4x₂ + 4x₃ + 0x₄ = 21 using Gaussian elimination with partial pivoting. Matrix: [0 4 2 8 | 24], [4 10 5 4 | 32], [4 5 6.5 2 | 26], [9 4 4 0 | 21]. Pivot a₁₁ = 0, so swap R₁ with R₄ (largest in column 1) → [9 4 4 0 | 21], [4 10 5 4 | 32], [4 5 6.5 2 | 26], [0 4 2 8 | 24]. Stage I: Divide R₁ by 9, eliminate x₁ from R₂, R₃, R₄. Continue elimination, choosing new pivots in each column sequentially (8.2222, 3.4594, 5.999). Final upper triangular system yields x₁ = 1.0, x₂ = 1.0, x₃ = 2.0, x₄ = 2.0.

📌 Example: Solve 3x + y - z = 3, 2x - 8y + z = -5, x - 2y + 9z = 8. Augmented matrix: [3 1 -1 | 3], [2 -8 1 | -5], [1 -2 9 | 8]. Stage I: R₂ ← R₂ - (2/3)R₁, R₃ ← R₃ - (1/3)R₁ → [3 1 -1 | 3], [0 -26/3 5/3 | -7], [0 -7/3 28/3 | 7]. Pivot in column 2 is -26/3 (non-zero). R₃ ← R₃ - (7/26)R₂ → [3 1 -1 | 3], [0 -26/3 5/3 | -7], [0 0 693/78 | 231/26]. Stage II: z = 1. y = 1. x = 1. Solution: x = y = z = 1.

📌 Example: Solve 3.15x - 1.96y + 3.85z = 12.95, 2.13x + 5.12y - 2.89z = -8.61, 5.92x + 3.05y + 2.15z = 6.88. Augmented matrix: [3.15 -1.96 3.85 | 12.95], [2.13 5.12 -2.89 | -8.61], [5.92 3.05 2.15 | 6.88]. R₂ ← R₂ - (2.13/3.15)R₁, R₃ ← R₃ - (5.92/3.15)R₁ → [3.15 -1.96 3.85 | 12.95], [0 6.4453 -5.4933 | -17.3667], [0 6.7335 -5.0855 | -17.4578]. Pivot = 6.4453. R₃ ← R₃ - (6.7335/6.4453)R₂ → [3.15 -1.96 3.85 | 12.95], [0 6.4453 -5.4933 | -17.3667], [0 0 0.6534 | 0.6853]. Back substitution: z = 0.6853/0.6534 = 1.0488215, y = (-17.3667 + 5.4933(1.0488215))/6.4453 = -1.8005692, x = (12.95 + 1.96(-1.8005692) - 3.85(1.0488215))/3.15 = 1.708864.

📌 Example: Solve x₁ + x₂ + x₃ + x₄ = 2, x₁ + x₂ + 3x₃ - 2x₄ = -6, 2x₁ + 3x₂ - x₃ + 2x₄ = 7, x₁ + 2x₂ + x₃ - x₄ = -2. Augmented matrix: [1 1 1 1 | 2], [1 1 3 -2 | -6], [2 3 -1 2 | 7], [1 2 1 -1 | -2]. Stage I: R₂ ← R₂ - R₁, R₃ ← R₃ - 2R₁, R₄ ← R₄ - R₁ → [1 1 1 1 | 2], [0 0 2 -3 | -8], [0 1 -3 0 | 3], [0 1 0 -2 | -4]. Pivot a₂₂ = 0, so swap R₂ and R₃ → [1 1 1 1 | 2], [0 1 -3 0 | 3], [0 0 2 -3 | -8], [0 1 0 -2 | -4]. R₄ ← R₄ - R₂ → [1 1 1 1 | 2], [0 1 -3 0 | 3], [0 0 2 -3 | -8], [0 0 3 -2 | -7]. Pivot = 2. R₄ ← R₄ - (3/2)R₃ → [1 1 1 1 | 2], [0 1 -3 0 | 3], [0 0 2 -3 | -8], [0 0 0 5/2 | 5]. Back substitution: x₄ = 2, x₃ = (-8 + 3(2))/2 = -1, x₂ = 3 + 3(-1) = 0, x₁ = 2 - 0 - (-1) - 2 = 1. Solution: x₁ = 1, x₂ = 0, x₃ = -1, x₄ = 2.

⭐ Key Takeaways

The Gaussian elimination method is a foundational direct method for solving linear systems, operating in two distinct phases: forward elimination to create an upper triangular matrix, followed by backward substitution. Pivoting, especially partial pivoting, is crucial to avoid zero pivots and mitigate round-off errors; it involves swapping rows to place the numerically largest element in the current column as the pivot. The method is applied systematically by first forming the augmented matrix, then performing elementary row operations with carefully computed multipliers. When a pivot is zero, the current row must be interchanged with a row below it that contains a non-zero (preferably the largest) element in that column. The solution is then extracted sequentially from the bottom up, substituting known values into each preceding equation.

🧠 Quick Revision Questions

1. What are the two main stages of the Gaussian elimination method?
2. Why is partial pivoting necessary, and what is the rule for choosing which row to swap with?
3. In the system x + y + z = 6, 2x + 3y + z = 10, 3x + 2y + 2z = 11, what is the first step (with pivoting) when using Gaussian elimination?
4. After forward elimination, how is the value of the last unknown found?
5. What distinguishes a direct method from an iterative method for solving linear systems?

📘 Lecture 11 — Solution of Linear System of Equations and Matrix Inversion Gauss–Jordon Elimination Method

📖 Overview: This lecture covers two powerful direct methods for solving systems of linear equations. The Gauss-Jordan elimination method extends Gaussian elimination to produce a diagonal matrix, directly yielding the solution without back substitution. Crout’s reduction method (also known as LU decomposition) factors the coefficient matrix into lower and upper triangular matrices, enabling efficient solution through forward and backward substitution. These techniques are fundamental for numerical analysis and computer programming.

🗂️ Topics Covered

The lecture begins with the Gauss-Jordan elimination method, which reduces a system to an equivalent diagonal form by making elements above and below the diagonal simultaneously zero. Several worked examples demonstrate the procedure for 3x3 and 4x4 systems with partial pivoting. The second half introduces Crout’s reduction method, which decomposes the coefficient matrix [A] into a lower-triangular matrix [L] and an upper-triangular matrix [U] with 1’s on the main diagonal. Detailed step-by-step examples show how to compute the [L] and [U] matrices, then solve for unknowns using forward substitution (Ly = B) and backward substitution (Ux = y).

📝 Lecture Summary

Gauss–Jordon Elimination Method

This method is a variation of Gaussian elimination. In this method, the elements above and below the diagonal are simultaneously made zero. That is, a given system is reduced to an equivalent diagonal form using elementary transformations. Then the solution of the resulting diagonal system is obtained. Sometimes, we normalize the pivot row with respect to the pivot element before elimination. Partial pivoting is also used whenever the pivot element becomes zero.

💡 Why this matters: Unlike Gaussian elimination, Gauss-Jordan eliminates the need for back substitution — the solution appears directly as the last column of the augmented matrix.

Example
Solve the system of equations using Gauss-Jordan elimination method:

x + 2y + z = 8
2x + 3y + 4z = 20
4x + 3y + 2z = 16

Solution
In matrix notation, the given system can be written as:

[1  2  1] [x]   [8 ]
[2  3  4] [y] = [20]
[4  3  2] [z]   [16]

After row operations (-2)R1+R2 and (-4)R1+R3:

[1  2  1] [x]   [8 ]
[0 -1  2] [y] = [4 ]
[0 -5 -2] [z]   [-16]

Eliminate y from the first and third rows using the second row:

[1  0  5] [x]   [16]
[0 -1  2] [y] = [4 ]
[0  0 -12][z]   [-36]

Normalize the third row (divide by -12):

[1  0  5] [x]   [16]
[0 -1  2] [y] = [4 ]
[0  0  1] [z]   [3 ]

Eliminate z from the first and second rows:

[1  0  0] [x]   [1 ]
[0 -1  0] [y] = [-2]
[0  0  1] [z]   [3 ]

The solution is x = 1, y = 2, z = 3.

Example
Solve using Gauss-Jordan method:

10x + y + z = 12
2x + 10y + z = 13
x + y + 5z = 7

The augmented matrix is:

[10  1  1 | 12]
[2  10  1 | 13]
[1   1  5 | 7 ]

After multiple row operations (R1 - 9R3, then R2 - 2R1 and R3 - R1, then R2 - 3R3, then R1 + 8R2 and R3 - 9R2, then dividing R3 by -473, then R1 - 420R3 and R2 - 58R3), the system reduces to reduced echelon form giving x = y = z = 1.

Example
Solve using Gauss-Jordan method:

10x₁ + x₂ + x₃ = 12
x₁ + 10x₂ - x₃ = 10
x₁ - 2x₂ + 10x₃ = 9

The solution gives x₁ = x₂ = x₃ = 1.

Example
Solve using Gauss-Jordan method:

x + 2y + z - w = -2
2x + 3y - z + 2w = 7
x + y + 3z - 2w = -6
x + y + z + w = 2

Through systematic row operations, the final reduced matrix gives:
x = 1, y = 0, z = -1, w = 2.

Crout’s Reduction Method

Here the coefficient matrix [A] of the system of equations is decomposed into the product of two matrices [L] and [U], where [L] is a lower-triangular matrix and [U] is an upper-triangular matrix with 1’s on its main diagonal.

For a general 3x3 matrix:

[l₁₁  0   0 ] [1  u₁₂  u₁₃]   [a₁₁  a₁₂  a₁₃]
[l₂₁  l₂₂ 0 ] [0  1   u₂₃] = [a₂₁  a₂₂  a₂₃]
[l₃₁  l₃₂  l₃₃] [0  0    1 ]   [a₃₁  a₃₂  a₃₃]

The sequence of steps:

1. Multiply all rows of [L] by the first column of [U]: l₁₁ = a₁₁, l₂₁ = a₂₁, l₃₁ = a₃₁ (first column of [L] = first column of [A]).
2. Multiply the first row of [L] by the second and third columns of [U]: u₁₂ = a₁₂/l₁₁, u₁₃ = a₁₃/l₁₁.
3. Multiply the 2nd and 3rd rows of [L] by the 2nd column of [U]: l₂₂ = a₂₂ - l₂₁u₁₂, l₃₂ = a₃₂ - l₃₁u₁₂.
4. Multiply the 2nd row of [L] by the 3rd column of [U]: u₂₃ = (a₂₃ - l₂₁u₁₃)/l₂₂.
5. Multiply the 3rd row of [L] by the 3rd column of [U]: l₃₃ = a₃₃ - l₃₁u₁₃ - l₃₂u₂₃.

This algorithm is generalized to any linear system of order n.

To solve A = (B): Let [A] = [L][U], then [L]U = (B). Substituting U = (Z), we get L = (B). Having computed Z (forward substitution), compute X from U = (Z) (back substitution).

This method is also known as the Cholesky reduction method. It is very popular from a computer programming point of view, since the storage space reserved for matrix [A] can be used to store the elements of [L] and [U] at the end of computation. This method fails if any aᵢᵢ = 0.

💡 Why this matters: LU decomposition is particularly efficient when solving multiple systems with the same coefficient matrix but different right-hand sides — the decomposition is done once.

Example
Solve using Crout’s reduction method:

5x₁ - 2x₂ + x₃ = 4
7x₁ + x₂ - 5x₃ = 8
3x₁ + 7x₂ + 4x₃ = 10

Solution
Let [L][U] = [A]:

[l₁₁ 0  0 ] [1  u₁₂  u₁₃]   [5  -2  1 ]
[l₂₁ l₂₂ 0 ] [0  1   u₂₃] = [7   1  -5]
[l₃₁ l₃₂ l₃₃] [0  0    1 ]   [3   7   4 ]

Step I: l₁₁ = 5, l₂₁ = 7, l₃₁ = 3
Step II: u₁₂ = -2/5, u₁₃ = 1/5
Step III: l₂₂ = 1 - 7(-2/5) = 19/5, l₃₂ = 7 - 3(-2/5) = 41/5
Step IV: u₂₃ = (-5 - 7(1/5))/(19/5) = -32/19
Step V: l₃₃ = 4 - 3(1/5) - (41/5)(-32/19) = 327/19

System becomes [L][U][X] = (B):

[5   0    0  ] [1  -2/5  1/5 ] [x₁]   [4 ]
[7  19/5  0  ] [0  1   -32/19] [x₂] = [8 ]
[3  41/5 327/19] [0  0     1  ] [x₃]   [10]

Let U = (Z), then L = (B):

5z₁ = 4 → z₁ = 4/5
7z₁ + (19/5)z₂ = 8 → z₂ = 12/19
3z₁ + (41/5)z₂ + (327/19)z₃ = 10 → z₃ = 46/327

Back substitution from U = (Z):

x₃ = 46/327
x₂ - (32/19)x₃ = 12/19 → x₂ = 284/327
x₁ - (2/5)x₂ + (1/5)x₃ = 4/5 → x₁ = 366/327

The solution is x₁ = 366/327, x₂ = 284/327, x₃ = 46/327.

Example
Solve using Crout’s method:

2x - 3y + 10z = 3
-x + 4y + 2z = 20
5x + 2y + z = -12

Here [L] has 1’s on the diagonal (different convention):

[1  0  0] [u₁₁  u₁₂  u₁₃]   [2  -3  10]
[l₂₁ 1  0] [0   u₂₂  u₂₃] = [-1  4   2 ]
[l₃₁ l₃₂ 1] [0   0   u₃₃]   [5   2   1 ]

Computing: u₁₁=2, u₁₂=-3, u₁₃=10, l₂₁=-1/2, l₃₁=5/2, u₂₂=5/2, u₂₃=7, l₃₂=19/5, u₃₃=-253/5

Let UX = Y, then LY = B:

y₁ = 3
(-1/2)y₁ + y₂ = 20 → y₂ = 43/2
(5/2)y₁ + (19/5)y₂ + y₃ = -12 → y₃ = -506/5

From UX = Y:

2x - 3y + 10z = 3
(5/2)y + 7z = 43/2
(-253/5)z = -506/5

By back substitution: x = -4, y = 3, z = 2.

Example
Solve using Crout’s method:

x + 3y + 8z = 4
x + 4y + 3z = -2
x + 3y + 4z = 1

Computing: u₁₁=1, u₁₂=3, u₁₃=8, l₂₁=1, l₃₁=1, u₂₂=1, u₂₃=-5, l₃₂=0, u₃₃=-4

From LY = B: y₁=4, y₂=-2, y₃=-3
From UX = Y: x+3y+8z=4, y-5z=-2, -4z=-3
Solution: z = 3/4, y = 7/4, x = 29/4.

⭐ Key Takeaways

The Gauss-Jordan elimination method reduces a system directly to a diagonal form by making elements both above and below the diagonal zero simultaneously, eliminating the need for back substitution — this is its key advantage over Gaussian elimination. Crout’s reduction method (LU decomposition) factors [A] into [L] (lower triangular) and [U] (upper triangular with 1’s on the diagonal), then solves through forward substitution ([L]Z = B) followed by back substitution ([U]X = Z), making it highly efficient for computer implementation where the same matrix is used with multiple right-hand sides. In Crout’s method, the first column of [L] equals the first column of [A], and the first row of [U] (except the first element which is 1) is obtained by dividing by l₁₁ — failure occurs if any diagonal element becomes zero.

🧠 Quick Revision Questions

1. What is the fundamental difference between Gauss-Jordan elimination and Gaussian elimination in terms of the final form of the matrix?
2. In Crout’s reduction method, what are the values on the main diagonal of the [U] matrix?
3. For the system represented by [L]U = (B), what two substitution steps are performed and in what order?
4. What condition causes Crout’s reduction method to fail?
5. In the Gauss-Jordan method for the system {x + 2y + z = 8, 2x + 3y + 4z = 20, 4x + 3y + 2z = 16}, what was the purpose of normalizing the third row (dividing by -12) before eliminating z from the first and second rows?

📘 Lecture 12 — Solution of Linear System of Equations and Matrix Inversion: Jacobi’s Method

📖 Overview: This lecture introduces Jacobi’s Method, an iterative technique for solving linear systems of equations, particularly advantageous when the coefficient matrix is sparse. It covers the method’s derivation, the crucial condition of diagonal dominance for convergence, and demonstrates the iterative process through multiple worked examples.

🗂️ Topics Covered

The lecture begins by introducing Jacobi’s Method as an iterative approach for solving linear systems, especially sparse ones. It then states the requirement for the coefficient matrix to be strictly diagonally dominant. The iterative formulas are derived from the original system by solving for each variable. Two detailed examples are presented, showing the step-by-step calculation of successive approximations from an initial guess of zero. A third example demonstrates how to rearrange a non-diagonally-dominant system to meet the convergence condition before applying the method.

📝 Lecture Summary

Solution of Linear System of Equations and Matrix Inversion

This section introduces Jacobi’s Method as an iterative method, where an initial approximate solution is assumed and progressively improved. This is in contrast to direct methods like Gaussian elimination.

💡 Why this matters: Iterative methods are highly efficient for sparse matrices (matrices with many zero elements), which are common in fields like solving partial differential equations. They require less computer memory than direct methods for such systems.

Jacobi’s Method

The method begins by assuming the coefficient matrix [A] is strictly diagonally dominant. This means that for every row, the absolute value of the diagonal element is greater than the sum of the absolute values of all other elements in that row.

|a_ii| > Σ |a_ij| for j ≠ i

It is also assumed that no diagonal element is zero. If it is, the equations should be rearranged.

The system of n equations is then rewritten to isolate each variable on the left-hand side:

x_i = (b_i - Σ a_ij * x_j) / a_ii for j ≠ i

Given an initial approximation vector (x₁⁽¹⁾, x₂⁽¹⁾, …, xₙ⁽¹⁾)ᵀ, the second approximation is computed by substituting the first into the right-hand side. This process yields the general formula for the (r+1)th approximation:

x_i^(r+1) = (b_i - Σ a_ij * x_j^(r)) / a_ii for j ≠ i

This is known as the method of simultaneous displacements because all variables from the r-th iteration are used to compute all new (r+1)-th values in one step, without updating them during the iteration.

🔑 Definition — Diagonal Dominance: A matrix condition where the diagonal element in each row is larger in magnitude than the sum of the magnitudes of all other elements in that row. It is a sufficient condition for the convergence of Jacobi’s method.

📐 Formula: x_i^(r+1) = (b_i - Σ a_ij * x_j^(r)) / a_ii for j ≠ i → The new value of variable x_i is calculated by dividing the constant b_i minus the sum of all off-diagonal terms (from row i, each multiplied by their current iteration r value) by the diagonal element a_ii.

Example 1

Problem: Solve the system using Jacobi’s iterative method for five iterations.

83x + 11y - 4z = 95 7x + 52y + 13z = 104 3x + 8y + 29z = 71

Solution: First, the equations are rewritten.

x = (95/83) - (11/83)y + (4/83)z → x = 1.1446 - 0.1325y + 0.0482z y = (104/52) - (7/52)x - (13/52)z → y = 2.0 - 0.1346x - 0.25z z = (71/29) - (3/29)x - (8/29)y → z = 2.4483 - 0.1035x - 0.2759y

The initial approximation is (x, y, z) = (0, 0, 0).

Iteration 1: Substitute the initial guess into the right-hand side.

x = 1.1446, y = 2.0, z = 2.4483 → (1.1446, 2.0, 2.4483)

Iteration 2: Substitute the values from Iteration 1.

x = 1.1446 - 0.1325(2.0) + 0.0482(2.4483) = 0.9976 y = 2.0 - 0.1346(1.1446) - 0.25(2.4483) = 1.2339 z = 2.4483 - 0.1035(1.1446) - 0.2759(2.0) = 1.7424 → (0.9976, 1.2339, 1.7424)

Further iterations yield the sequence:

The values appear to be converging towards the exact solution.

Example 2

Problem: Solve the system using Jacobi’s method (four iterations).

8x - 3y + 2z = 20 4x + 11y - z = 33 6x + 3y + 12z = 35

Solution: The system is diagonally dominant. The equations are rewritten.

x = (1/8)[20 + 3y - 2z] y = (1/11)[33 - 4x + z] z = (1/12)[35 - 6x - 3y]

Starting from (0, 0, 0):

Iteration 1: x₁ = 2.5, y₁ = 3, z₁ = 2.916667 Iteration 2: x₂ = 2.895833, y₂ = 2.3560606, z₂ = 0.9166666 Iteration 3: x₃ = 3.1543561, y₃ = 2.030303, z₃ = 0.8797348 Iteration 4: x₄ = 3.0419299, y₄ = 1.9329373, z₄ = 0.8319128

Example 3

Problem: Solve the system using Jacobi’s method (four iterations).

3x + 4y + 15z = 54.8 x + 12y + 3z = 39.66 10x + y - 2z = 7.74

Solution: The original system is not diagonally dominant. It is rearranged to meet the condition.

10x + y - 2z = 7.74 x + 12y + 3z = 39.66 3x + 4y + 15z = 54.8

Now, the system is diagonally dominant. The equations are rewritten.

x = (1/10)[7.74 - y + 2z] y = (1/12)[39.66 - x - 3z] z = (1/15)[54.8 - 3x - 4y]

Starting from (0, 0, 0): Iteration 1: x₁ = 0.774, y₁ = 3.305, z₁ = 3.6533333 Iteration 2: x₂ = 1.3908333, y₂ = 2.3271667, z₂ = 3.1949778 Iteration 3: x₃ = 1.1802789, y₃ = 2.3903528, z₃ = 2.7545889 Iteration 4: x₄ = 1.0781704, y₄ = 2.51779962, z₄ = 2.7798501

⭐ Key Takeaways

1. Jacobi’s method is an iterative technique that solves linear systems by repeatedly substituting the entire previous approximation into a set of reformulated equations, making it a method of simultaneous displacements.
2. A sufficient condition for the convergence of the Jacobi method is that the coefficient matrix is strictly diagonally dominant. If a system is not diagonally dominant, it must be rearranged before the method can be reliably applied.
3. The method is particularly advantageous for sparse matrices (where many elements are zero), as it is more memory-efficient than direct methods like Gaussian elimination.
4. The process begins with an initial guess (often the zero vector) and generates a sequence of approximations that, if convergent, approach the true solution to the system.

🧠 Quick Revision Questions

1. What is the key mathematical condition that is sufficient for the Jacobi iterative method to converge?
2. Explain why Jacobi’s method is called the "method of simultaneous displacements."
3. Given the system 10x - y + z = 12 and 2x + 10y - z = 11 and -x + y + 10z = 10, write the equations in the form suitable for Jacobi's method.
4. For the system in question 3, perform the first iteration starting from (0,0,0). What is the resulting approximate solution vector?
5. If a system of linear equations is not diagonally dominant, what must be done before applying Jacobi’s method, and why?

📘 Lecture 13 — Solution of Linear System of Equations and Matrix Inversion

📖 Overview: This lecture introduces the Gauss-Seidel Iteration Method for solving systems of linear equations. It explains how this method differs from Jacobi's method by using updated values immediately, making it faster and more efficient for large systems.

🗂️ Topics Covered

The lecture covers the Gauss-Seidel iteration method formulation, its comparison with Jacobi's method, the general compact form of the procedure, and multiple solved examples demonstrating the iterative process for various systems of equations including 3×3 and 4×4 systems.

📝 Lecture Summary

Solution of Linear System of Equations and Matrix Inversion – Gauss–Seidel Iteration Method

The Gauss-Seidel Iteration Method is an iterative technique for solving linear systems of the form a₁₁x₁ + a₁₂x₂ + … + a₁ₙxₙ = b₁, a₂₁x₁ + a₂₂x₂ + … + a₂ₙxₙ = b₂, … , aₙ₁x₁ + aₙ₂x₂ + … + aₙₙxₙ = bₙ.

In Jacobi’s method, the (r+1)th approximation uses only values from the rth iteration. No element of the new approximation replaces the old one during computation. In contrast, in the Gauss-Seidel method, the corresponding elements of x^(r+1) replace those of x^(r) as soon as they become available. Hence, it is called the method of successive displacements.

🔑 Definition — Gauss-Seidel Iteration: For a system of n linear equations, the (r+1)th approximation is computed using the most recently updated values. For the first equation, we use x₂^(r), x₃^(r), …, xₙ^(r). For the second equation, we use x₁^(r+1), x₃^(r), …, xₙ^(r). For the third, we use x₁^(r+1), x₂^(r+1), x₄^(r), …, and so on.

📐 Formula: For each variable xᵢ at iteration (r+1): xᵢ^(r+1) = (bᵢ - Σⱼ₌₁^(i-1) aᵢⱼ xⱼ^(r+1) - Σⱼ₌ᵢ₊₁ⁿ aᵢⱼ xⱼ^(r)) / aᵢᵢ

This means: To find each new value, subtract the contributions of all previously updated variables and the remaining old variables from the constant term, then divide by the diagonal coefficient.

💡 Why this matters: The immediate replacement of values in Gauss-Seidel typically leads to faster convergence than Jacobi's method, requiring fewer iterations to reach the same accuracy.

Example — System of 4 Equations

Solve the system using Gauss-Seidel method with first five iterations: 4x₁ - x₂ - x₃ = 2 -x₁ + 4x₂ - x₄ = 2 -x₁ + 4x₃ - x₄ = 1 -x₂ - x₃ + 4x₄ = 1

📌 Example: The system is rewritten as: x₁ = 0.5 + 0.25x₂ + 0.25x₃ x₂ = 0.5 + 0.25x₁ + 0.25x₄ x₃ = 0.25 + 0.25x₁ + 0.25x₄ x₄ = 0.25 + 0.25x₂ + 0.25x₃

Starting with x₂ = x₃ = x₄ = 0 in the first equation: x₁^(1) = 0.5 Using x₁^(1) = 0.5 and x₃ = x₄ = 0: x₂^(1) = 0.5 + (0.25)(0.5) + 0 = 0.625 Using x₁^(1) = 0.5, x₂^(1) = 0.625, and x₄ = 0: x₃^(1) = 0.25 + (0.25)(0.5) + 0 = 0.375 Using x₂^(1) = 0.625 and x₃^(1) = 0.375: x₄^(1) = 0.25 + (0.25)(0.625) + (0.25)(0.375) = 0.5

Continuing iterations: Iteration 1: x₁=0.5, x₂=0.625, x₃=0.375, x₄=0.5 Iteration 2: x₁=0.75, x₂=0.8125, x₃=0.5625, x₄=0.59375 Iteration 3: x₁=0.84375, x₂=0.85938, x₃=0.60938, x₄=0.61719 Iteration 4: x₁=0.86719, x₂=0.87110, x₃=0.62110, x₄=0.62305 Iteration 5: x₁=0.87305, x₂=0.87402, x₃=0.62402, x₄=0.62451

Example — 3×3 System with Diagonal Dominance

Solve by Gauss-Seidel iterative method (four iterations): 8x - 3y + 2z = 20 4x + 11y - z = 33 6x + 3y + 12z = 35

📌 Example: The system is diagonally dominant (each diagonal coefficient is greater than the sum of absolute values of off-diagonal coefficients in its row). Rewrite: x = (1/8)[20 + 3y - 2z] y = (1/11)[33 - 4x + z] z = (1/12)[35 - 6x - 3y]

Starting with x₀ = y₀ = z₀ = 0: Iteration 1: x₁ = (1/8)[20] = 2.5, y₁ = (1/11)[33 - 4(2.5) + 0] = 2.0909091, z₁ = (1/12)[35 - 6(2.5) - 3(2.0909091)] = 1.1439394 Iteration 2: x₂ = (1/8)[20 + 3(2.0909091) - 2(1.1439394)] = 2.9981061, y₂ = (1/11)[33 - 4(2.9981061) + 1.1439394] = 2.0137741, z₂ = (1/12)[35 - 6(2.9981061) - 3(2.0137741)] = 0.9141701 Iteration 3: x₃ = 3.0266228, y₃ = 1.9825163, z₃ = 0.9077262 Iteration 4: x₄ = 3.0165121, y₄ = 1.9856071, z₄ = 0.8319128

Example — System Requiring Equation Reordering

Solve using Gauss-Seidel iteration method: 28x + 4y - z = 32 x + 3y + 10z = 24 2x + 17y + 4z = 35

📌 Example: The original system is not diagonally dominant. The equations must be rearranged to make the system diagonally dominant: 28x + 4y - z = 32 (keep as is) 2x + 17y + 4z = 35 (move to second position) x + 3y + 10z = 24 (move to third position)

Now rewrite: x = (1/28)[32 - 4y + z] y = (1/17)[35 - 2x - 4z] z = (1/10)[24 - x - 3y]

Starting with y = z = 0: x₁ = 32/28 = 1.1428571 Using x₁ = 1.1428571, z = 0: y₁ = (1/17)[35 - 2(1.1428571)] = 1.9243697 Using x₁ = 1.1428571, y₁ = 1.9243697: z₁ = (1/10)[24 - 1.1428571 - 3(1.9243697)] = 1.7084034

Continuing iterations: Iteration 2: x₂ = 0.9289615, y₂ = 1.5475567, z₂ = 1.8408368 Iteration 3: x₃ = 0.9875932, y₃ = 1.5090274, z₃ = 1.8485325 Iteration 4: x₄ = 0.9933008, y₄ = 1.5070158, z₄ = 1.8485652

Example — 4×4 System

Solve using Gauss-Seidel iteration method (four iterations): 10x - 2y - z - w = 3 -2x + 10y - z - w = 15 -x - y + 10z - 2w = 27 -x - y - 2z + 10w = -9

📌 Example: The system is diagonally dominant. Rewrite: x = (1/10)[3 + 2y + z + w] y = (1/10)[15 + 2x + z + w] z = (1/10)[27 + x + y + 2w] w = (1/10)[-9 + x + y + 2z]

Starting with y = z = w = 0: x₁ = 0.3 Using x₁ = 0.3: y₁ = (1/10)[15 + 2(0.3)] = 1.56 Using x₁ = 0.3, y₁ = 1.56: z₁ = (1/10)[27 + 0.3 + 1.56] = 2.886 Using x₁ = 0.3, y₁ = 1.56, z₁ = 2.886: w₁ = (1/10)[-9 + 0.3 + 1.56 + 2(2.886)] = -0.1368

Continuing iterations: Iteration 2: x₂ = 0.88692, y₂ = 1.952304, z₂ = 2.9565624, w₂ = -0.0247651 Iteration 3: x₃ = 0.9836405, y₃ = 1.9899087, z₃ = 2.9924019, w₃ = -0.0041647 Iteration 4: x₄ = 0.9968054, y₄ = 1.9981848, z₄ = 2.9986661, w₄ = -0.0007677

When to Stop the Iterative Process

We stop the iterative process when we achieve the required accuracy. For example, if asked to find values accurate to six decimal places, and consecutive values are 1.919326355, 1.919327145, 1.919327128, the accuracy up to seven decimal places is achieved, so we stop. The stopping criterion depends on the desired precision specified in the problem.

⭐ Key Takeaways

The Gauss-Seidel iteration method is an improvement over Jacobi's method because it uses the most recently updated values immediately, leading to faster convergence. The system must be diagonally dominant for convergence, and if it isn't, the equations must be rearranged to achieve diagonal dominance. The general formula for each variable uses the sum of products of coefficients with previously updated values and remaining old values. The iterative process continues until the desired accuracy is reached, which is determined by comparing consecutive approximations. Gauss-Seidel is particularly effective for large sparse systems where direct methods are computationally expensive.

🧠 Quick Revision Questions

1. What is the key difference between Jacobi's method and Gauss-Seidel method in terms when updated values are used?
2. What condition must a system of equations satisfy for the Gauss-Seidel method to converge?
3. For a 4×4 system using Gauss-Seidel, which variable values are used to compute x₂^(r+1)?
4. In the example 10x - 2y - z - w = 3, what is the iterative formula for z?
5. How do you determine when to stop the Gauss-Seidel iterative process?

📘 Lecture 14 — Solution of Linear System of Equations and Matrix Inversion Relaxation Method

📖 Overview: This lecture introduces the Relaxation Method, an iterative technique for solving linear systems of equations. It is based on systematically reducing the largest residual (error) in each equation to zero by updating the corresponding variable. This method is attributed to Southwell and is particularly useful for achieving fast convergence when equations are properly ordered.

🗂️ Topics Covered

The lecture begins by defining the residual of each equation in a linear system and establishing the core principle of the relaxation method: iteratively reducing the largest residual to zero. It covers the procedure for reordering equations so that the largest diagonal coefficients appear on the diagonal. The lecture then details the formula for calculating the increment to a variable based on its residual and diagonal coefficient. Finally, several worked examples are presented to demonstrate the step-by-step iterative process, from initial residuals to final solution convergence.

📝 Lecture Summary

Relaxation Method

This is an iterative method attributed to Southwell. Consider the system of equations: [ a_{11}x_1 + a_{12}x_2 + \dots + a_{1n}x_n = b_1 ] [ a_{21}x_1 + a_{22}x_2 + \dots + a_{2n}x_n = b_2 ] [ \vdots ] [ a_{n1}x_1 + a_{n2}x_2 + \dots + a_{nn}x_n = b_n ]

Let ( X^{(p)} = (x_1^{(p)}, x_2^{(p)}, \dots, x_n^{(p)})^T ) be the solution vector after the ( p )-th iteration.

If ( R_i^{(p)} ) denotes the residual of the ( i )-th equation, defined by: [ R_i^{(p)} = b_i - a_{i1}x_1^{(p)} - a_{i2}x_2^{(p)} - \dots - a_{in}x_n^{(p)} ] We can improve the solution vector successively by reducing the largest residual to zero at that iteration. This is the basic idea of the relaxation method.

To achieve fast convergence, we take all terms to one side and then reorder the equations so that the largest negative coefficients in the equations appear on the diagonal.

Now, if at any iteration, ( R_i ) is the largest residual in magnitude, then we give an increment to ( x_i ); ( a_{ii} ) being the coefficient of ( x_i ): [ dx_i = \frac{R_i}{a_{ii}} ] In other words, we change ( x_i ) to ( (x_i + dx_i) ) to relax ( R_i ), i.e., to reduce ( R_i ) to zero.

🔑 Definition — Residual: The residual ( R_i ) of the ( i )-th equation is the difference between the constant term ( b_i ) and the current computed value of the left-hand side. It measures the error in that equation at a given iteration. 📐 Formula: ( R_i = b_i - \sum_{j=1}^{n} a_{ij}x_j ) and ( dx_i = R_i / a_{ii} ) → The increment to a variable equals its residual divided by its diagonal coefficient.

Example 1

Solve the system of equations: [ 6x_1 - 3x_2 + x_3 = 11 ] [ 2x_1 + x_2 - 8x_3 = -15 ] [ x_1 - 7x_2 + x_3 = 10 ] by the relaxation method, starting with the vector (0, 0, 0).

Solution: At first, transfer all terms to the right-hand side and reorder the equations so that the largest coefficients appear on the diagonal: [ 0 = 11 - 6x_1 + 3x_2 - x_3 ] [ 0 = 10 - x_1 + 7x_2 - x_3 ] [ 0 = -15 - 2x_1 - x_2 + 8x_3 ] after interchanging the 2nd and 3rd equations.

Starting with the initial solution vector (0, 0, 0), we find the residuals: [ R_1 = 11, R_2 = 10, R_3 = -15 ] The largest residual in magnitude is ( R_3 = -15 ). We introduce a change, ( dx_3 ) in ( x_3 ): [ dx_3 = \frac{R_3}{a_{33}} = \frac{15}{8} = 1.875 ]

The iterative process continues until all residuals are very small. The iteration table is:

The numerical solution is ( x_1 = 1.0017, x_2 = -0.9901, x_3 = 2.0017 ). The exact solution is ( x_1 = 1.0, x_2 = -1.0, x_3 = 2.0 ).

📌 Example: In iteration 0, the largest residual is -15 for equation 3. The increment for ( x_3 ) is calculated as ( dx_3 = -(-15)/8 = 15/8 = 1.875 ). After updating ( x_3 ) to 1.875, new residuals are computed for iteration 1: ( R_1 = 11 - 6(0) + 3(0) - 1(1.875) = 9.125 ), ( R_2 = 10 - 1(0) + 7(0) - 1(1.875) = 8.125 ), ( R_3 = -15 - 2(0) - 1(0) + 8(1.875) = 0 ). The process continues until all residuals are small.

Example 2

Solve by relaxation method, the equation: [ 10x - 2y - 2z = 6 ] [ -x - 10y - 2z = 7 ] [ -x - y + 10z = 8 ]

Solution: The residuals ( r_1, r_2, r_3 ) are given by: [ r_1 = 6 - 10x + 2y + 2z ] [ r_2 = 7 + x - 10y + 2z ] [ r_3 = 8 + x + y - 10z ]

The operation table is:

The relaxation table is:

Explanation: (1) In L4, the largest residual is 8. To reduce it, we give an increment of ( 8/10 = 0.8 \approx 1 ) to ( z ). The resulting residuals are obtained by L4 + (1)L3, i.e., line L5. (2) In line L5, the largest residual is 9. Increment = ( 9/10 = 0.9 \approx 1 ) to ( y ). The resulting residuals (= L6) are L5 + 1.L2. (3) In line L6, the largest residual is 10. Increment = ( 10/10 \approx 1 ) to ( x ). The resulting residuals (= L7) are L6 + 1.L1.

The exact solution is ( x = 1, y = 1, z = 1 ).

Example 3

Solve the system by relaxation method: [ 9x - y + 2z = 7 ] [ x + 10y - 2z = 15 ] [ 2x - 2y - 13z = -17 ]

Solution: The residuals ( r_1, r_2, r_3 ) are given by: [ r_1 = 9 - 9x + y - 2z ] [ r_2 = 15 - x - 10y + 2z ] [ r_3 = -17 - 2x + 2y + 13z ]

Operation table:

Relaxation table:

Then: [ x = 0.89 + 0.028 = 0.918 ] [ y = 1 + 0.61 + 0.039 = 1.649 ] [ z = 1 + 0.19 + 0.00523 = 1.19523 ]

Substituting the values: [ r_1 = 9 - 9(0.918) + 1.649 - 2(1.19523) = -0.00346 ] [ r_2 = 15 - 0.918 - 10(1.649) + 2(1.19523) = -0.1754 ] [ r_3 = -17 - 2(0.918) + 2(1.649) + 13(1.19523) = -0.00001 ] Which agrees with the final residuals.

💡 Why this matters: The relaxation method provides a systematic way to iteratively solve linear systems by focusing on the equation with the largest error, making it efficient for certain types of problems.

⭐ Key Takeaways

The relaxation method is an iterative technique that works by computing residuals for each equation and then systematically reducing the largest residual to zero by updating the corresponding variable. The increment for a variable is calculated as its residual divided by its diagonal coefficient. To improve convergence, equations should be reordered so that the largest coefficients appear on the diagonal. The process continues until all residuals become sufficiently small, indicating that the solution has been reached. This method is particularly effective when the system has a strong diagonal dominance.

🧠 Quick Revision Questions

1. What is a residual in the context of the relaxation method, and how is it calculated?
2. What is the formula for calculating the increment (dx_i) to a variable x_i during the relaxation process?
3. Why is it important to reorder equations in the relaxation method?
4. In Example 1, why was the first increment applied to x₃ instead of x₁ or x₂?
5. When can the iterative process in the relaxation method be stopped?

📘 Lecture 15 — Solution of Linear System of Equations and Matrix Inversion

📖 Overview: This lecture focuses on methods for finding the inverse of a matrix, which is essential for solving linear systems of equations in the form A = (B). It covers the Gaussian Elimination Method and the Gauss-Jordan Method, two fundamental techniques for computing matrix inverses through row operations and augmented matrices. Understanding these methods is critical for numerical analysis and solving real-world engineering problems.

🗂️ Topics Covered

This lecture introduces matrix inversion as a method for solving linear systems of equations. It explains the Gaussian Elimination Method in two stages: first reducing the matrix to an upper triangular form, then reducing it to an identity matrix to obtain the inverse. It also covers the Gauss-Jordan Method, which directly reduces the matrix to an identity matrix using elementary row operations. Several detailed examples are provided for both methods, including the importance of partial pivoting for accuracy.

📝 Lecture Summary

Solution of Linear System of Equations and Matrix Inversion

Consider a system of equations in the form A = (B). One way of writing its solution is in the form (X) = [A]⁻¹(B). Thus, the solution to the system can also be obtained if the inverse of the coefficient matrix [A] is known. The product of two square matrices is an identity matrix [A][B] = [I], then [B] = [A]⁻¹ and [A] = [B]⁻¹. Every square non-singular matrix will have an inverse. Gauss elimination and Gauss-Jordan methods are popular among many methods available for finding the inverse of a matrix.

Gaussian Elimination Method

In this method, if A is a given matrix, we first place an identity matrix, whose order is same as that of A, adjacent to A, forming an augmented matrix. The inverse of A is computed in two stages. In the first stage, A is converted into an upper triangular form, using the Gaussian elimination method. In the second stage, the upper triangular matrix is reduced to an identity matrix by row transformations. All these operations are also performed on the adjacently placed identity matrix. When A is transformed into an identity matrix, the adjacent matrix gives the inverse of A. To increase accuracy, it is essential to employ partial pivoting.

Example: Use the Gaussian elimination method to find the inverse of the matrix A = [[1,1,1],[4,3,-1],[3,5,3]].

Solution: First, place an identity matrix adjacent to the given matrix: [[1,1,1,1,0,0],[4,3,-1,0,1,0],[3,5,3,0,0,1]].

Stage I (Reduction to upper triangular form): In the 1st column, 4 is the largest element, so interchange R1 and R2 (partial pivoting): [[4,3,-1,0,1,0],[1,1,1,1,0,0],[3,5,3,0,0,1]].

Divide R1 by 4: [[1,3/4,-1/4,0,1/4,0],[1,1,1,1,0,0],[3,5,3,0,0,1]].

Perform R₂ – R₁ → R₂: [[1,3/4,-1/4,0,1/4,0],[0,1/4,5/4,1,-1/4,0],[3,5,3,0,0,1]].

Perform R₃ – 3R₁ → R₃: [[1,3/4,-1/4,0,1/4,0],[0,1/4,5/4,1,-1/4,0],[0,11/4,15/4,0,-3/4,1]].

For the second column pivot, max (1/4, 11/4) is 11/4, so interchange R2 and R3: [[1,3/4,-1/4,0,1/4,0],[0,11/4,15/4,0,-3/4,1],[0,1/4,5/4,1,-1/4,0]].

Divide R2 by the pivot a₂₂ = 11/4: [[1,3/4,-1/4,0,1/4,0],[0,1,15/11,0,-3/11,4/11],[0,1/4,5/4,1,-1/4,0]].

Perform R₃ – (1/4)R₂ → R₃: [[1,3/4,-1/4,0,1/4,0],[0,1,15/11,0,-3/11,4/11],[0,0,10/11,1,-2/11,-1/11]].

Finally, divide R3 by (10/11): [[1,3/4,-1/4,0,1/4,0],[0,1,15/11,0,-3/11,4/11],[0,0,1,11/10,-1/5,-1/10]]. This is the upper triangular form.

Stage II (Reduction to an identity matrix): Perform (-1/4)R₃ + R₁ and (-15/11)R₃ + R₂: [[1,3/4,0,11/40,1/5,-1/40],[0,1,0,-3/2,0,1/2],[0,0,1,11/10,-1/5,-1/10]].

Finally, perform R₁ – (3/4)R₂ → R₁: [[1,0,0,7/5,1/5,-2/5],[0,1,0,-3/2,0,1/2],[0,0,1,11/10,-1/5,-1/10]].

Thus, A⁻¹ = [[7/5, 1/5, -2/5], [-3/2, 0, 1/2], [11/10, -1/5, -1/10]].

Example: Find the inverse of [[2,1,1],[3,2,3],[1,4,9]] using Gauss elimination.

Solution: Augmented matrix: [[2,1,1,1,0,0],[3,2,3,0,1,0],[1,4,9,0,0,1]].

After row operations R₂ – (3/2)R₁ and R₃ – (1/2)R₁: [[2,1,1,1,0,0],[0,1/2,3/2,-3/2,1,0],[0,7/2,17/2,-1/2,0,1]].

Perform R₃ – 7R₂: [[2,1,1,1,0,0],[0,1/2,3/2,-3/2,1,0],[0,0,-2,10,-7,1]].

If the inverse matrix is [[x₁₁,x₁₂,x₁₃],[x₂₁,x₂₂,x₂₃],[x₃₁,x₃₂,x₃₃]], the system is equivalent to three systems. Solving by back substitution:

From the first system: 2x₁₁ + x₂₁ + x₃₁ = 1, (1/2)x₂₁ + (3/2)x₃₁ = -3/2, -2x₃₁ = 10 → x₃₁ = -5, then x₂₁ = 12, x₁₁ = -3.

From the second system: 2x₁₂ + x₂₂ + x₃₂ = 0, (1/2)x₂₂ + (3/2)x₃₂ = 1, -2x₃₂ = -7 → x₃₂ = 7/2, then x₂₂ = -17/2, x₁₂ = 5/2.

From the third system: 2x₁₃ + x₂₃ + x₃₃ = 0, (1/2)x₂₃ + (3/2)x₃₃ = 0, -2x₃₃ = 1 → x₃₃ = -1/2, then x₂₃ = 3/2, x₁₃ = -1/2.

The inverse matrix is (1/2)[[-6,5,-1],[24,-17,3],[-10,7,-1]].

Example: Find the inverse of [[4,1,2],[2,3,-1],[1,-2,2]] using Gauss elimination.

Solution: Augmented matrix: [[4,1,2,1,0,0],[2,3,-1,0,1,0],[1,-2,2,0,0,1]].

After row operations R₂ – (1/2)R₁ and R₃ – (1/4)R₁: [[4,1,2,1,0,0],[0,5/2,-2,-1/2,1,0],[0,-9/4,3/2,-1/4,0,1]].

This gives three systems to solve. Solving by back substitution yields the inverse matrix. The final inverse is (1/3)[[-4,6,7],[5,-8,-2],[7,-6,-10]].

Gauss - Jordan Method

This method is similar to the Gaussian elimination method, with the essential difference that stage I of reducing the given matrix to an upper triangular form is not needed. The given matrix can be directly reduced to an identity matrix using elementary row operations.

Example: Find the inverse of A = [[1,1,1],[4,3,-1],[3,5,3]] by Gauss-Jordan method.

Solution: Augmented matrix: [[1,1,1,1,0,0],[4,3,-1,0,1,0],[3,5,3,0,0,1]].

Perform R₂ – 4R₁ → R₂: [[1,1,1,1,0,0],[0,-1,-5,-4,1,0],[3,5,3,0,0,1]].

Perform R₃ – 3R₁ → R₃: [[1,1,1,1,0,0],[0,-1,-5,-4,1,0],[0,2,0,-3,0,1]].

Perform R₂ + R₁ → R₁ and R₃ + 2R₂ → R₃: [[1,0,-4,-3,1,0],[0,-1,-5,-4,1,0],[0,0,-10,-11,2,1]].

Divide the third row by -10: [[1,0,-4,-3,1,0],[0,-1,-5,-4,1,0],[0,0,1,11/10,-1/5,-1/10]].

Perform R₁ + 4R₃ → R₁ and R₂ + 5R₃ → R₂: [[1,0,0,7/5,1/5,-2/5],[0,-1,0,-3/2,0,1/2],[0,0,1,11/10,-1/5,-1/10]].

Multiply R2 by -1: [[1,0,0,7/5,1/5,-2/5],[0,1,0,3/2,0,-1/2],[0,0,1,11/10,-1/5,-1/10]]. Note: there's a sign error in the text; the correct inverse is A⁻¹ = [[7/5, 1/5, -2/5], [3/2, 0, -1/2], [11/10, -1/5, -1/10]].

💡 Why this matters: Both methods produce the same inverse matrix, but Gauss-Jordan is more direct as it eliminates the need for a separate back-substitution stage, though it requires more row operations.

⭐ Key Takeaways

The most critical concepts from this lecture are: matrix inversion via the augmented matrix method, where an identity matrix is placed alongside the given matrix and row operations transform the given matrix into an identity matrix. The Gaussian elimination method proceeds in two distinct stages: first reducing to upper triangular form, then to identity form. The Gauss-Jordan method accomplishes this in one continuous process. Partial pivoting (reordering rows to bring the largest element to the pivot position) is essential for numerical accuracy. The final inverse matrix appears in the columns originally occupied by the identity matrix.

🧠 Quick Revision Questions

1. What is the augmented matrix in the Gaussian elimination method for finding matrix inverse?
2. Why is partial pivoting important in matrix inversion methods?
3. In Gaussian elimination, what is the goal of Stage I and Stage II?
4. What is the key difference between the Gaussian elimination method and the Gauss-Jordan method?
5. Given an upper triangular augmented matrix from Stage I, how do you solve for the inverse matrix elements?

📘 Lecture 16 — Eigen Value Problems

📖 Overview: This lecture introduces eigenvalue problems for square matrices and presents numerical methods for computing eigenvalues and eigenvectors. It focuses on the Power Method for finding the largest eigenvalue and its corresponding eigenvector, which is essential when solving the characteristic polynomial becomes computationally prohibitive.

🗂️ Topics Covered

The lecture begins with the definition of eigenvalues and eigenvectors, then explains why numerical methods are preferred over solving the characteristic polynomial. It introduces the Power Method procedure for finding the largest eigenvalue, demonstrates the method through multiple examples with full iterations, and discusses the mathematical theory behind why the method converges. Finally, it mentions using the inverse matrix to find the smallest eigenvalue.

📝 Lecture Summary

Eigen Value Problems

Let [A] be an n x n square matrix. Suppose there exists a scalar λ and a vector X = (x₁, x₂, ..., xₙ)ᵀ such that A = λ(X). Then λ is the eigenvalue and X is the corresponding eigenvector of matrix [A]. This can also be written as A - λI = (O), representing a set of n homogeneous equations with non-trivial solution provided |A - λI| = 0. This determinant, on expansion, gives an n-th degree polynomial called the characteristic polynomial of [A], which has n roots.

Finding the roots of the characteristic equation is laborious, so we look for better computational methods. The lecture discusses Power Method and Jacobi's Method for real and real-symmetric matrices.

🔑 Definition — Eigenvalue: A scalar λ such that A = λ(X) for some non-zero vector X 🔑 Definition — Eigenvector: The vector X satisfying A = λ(X) for the corresponding eigenvalue λ 🔑 Definition — Characteristic Polynomial: The n-th degree polynomial obtained from expanding |A - λI| = 0

Power Method

To compute the largest eigenvalue and corresponding eigenvector of the system A = λ(X), where [A] is real, symmetric, or unsymmetric, the Power Method is widely used.

Procedure: Step 1: Choose the initial vector such that the largest element is unity. Step 2: The normalized vector v⁽⁰⁾ is pre-multiplied by matrix [A]. Step 3: The resultant vector is again normalized. Step 4: This iteration process continues, repeatedly pre-multiplying the new normalized vector by [A] until required accuracy is obtained.

At the end, u⁽ᵏ⁾ = [A]v⁽ᵏ⁻¹⁾ = qₖv⁽ᵏ⁾, where qₖ is the desired largest eigenvalue and v⁽ᵏ⁾ is the corresponding eigenvector.

📌 Example: Find the eigenvalue of largest modulus and associated eigenvector of matrix A = [[2,3,2],[4,3,5],[3,2,9]] using Power Method.

Starting with v⁽⁰⁾ = (1,1,1)ᵀ: 1st iteration: u⁽¹⁾ = [A]v⁽⁰⁾ = (7,12,14)ᵀ, normalize using 14: v⁽¹⁾ = (1/2, 6/7, 1)ᵀ

2nd iteration: u⁽²⁾ = (39/7, 67/7, 171/7)ᵀ ≈ (5.571, 9.571, 24.429)ᵀ, normalize by 12.2143: v⁽²⁾ ≈ (0.456140, 0.783626, 1.0)ᵀ

3rd iteration: u⁽³⁾ ≈ (5.263158, 9.175438, 11.935672)ᵀ, normalize by 11.935672: v⁽³⁾ ≈ (0.44096, 0.776874, 1.0)ᵀ

4th iteration: u⁽⁴⁾ ≈ (5.18814, 9.07006, 11.86036)ᵀ, normalize by 11.8636: v⁽⁴⁾ ≈ (0.437435, 0.764737, 1.0)ᵀ

5th iteration: u⁽⁵⁾ ≈ (5.16908, 9.04395, 11.84178)ᵀ, normalize by 11.84178: v⁽⁵⁾ ≈ (0.436512, 0.763732, 1.0)ᵀ

After rounding, the largest eigenvalue λ ≈ 11.84 and eigenvector X ≈ (0.44, 0.76, 1.00)ᵀ

📌 Example: Find first three iterations for matrix A = [[7,6,-3],[-12,-20,24],[-6,-12,16]]

With v⁽⁰⁾ = (1,1,1)ᵀ: 1st iteration: u⁽¹⁾ = (10,-8,-2)ᵀ, normalize by 10: v⁽¹⁾ = (1, -0.8, -0.2)ᵀ, q₁ = 10

2nd iteration: u⁽²⁾ = (2.8,-0.8,0.4)ᵀ, normalize by 2.8: v⁽²⁾ = (1, -0.2857, 0.1428)ᵀ, q₂ = 2.8

3rd iteration: u⁽³⁾ = (4.8574,-2.8588,-0.2868)ᵀ, normalize by 4.8574: v⁽³⁾ = (1, -0.5885, -0.0590)ᵀ

📌 Example: Find first three iterations for matrix A = [[1,-1,0],[-2,4,-2],[0,-1,2]] with x⁽⁰⁾ = (-1,2,1)ᵀ

1st iteration: u⁽¹⁾ = (-3,8,0)ᵀ, normalize by 8: x⁽¹⁾ = (-3/8, 1, 0)ᵀ, q₁ = 8

2nd iteration: u⁽²⁾ = (-1.375, 4.75, -1)ᵀ, normalize by 4.75: x⁽²⁾ = (-0.28947, 1, -0.2152)ᵀ, q₂ = 4.75

3rd iteration: u⁽³⁾ = (-1.28947, 4.99998, -1.42104)ᵀ, normalize by 4.99998: x⁽³⁾ = (-0.25789, 1, -0.28420)ᵀ

Mathematical Foundation

Let λ₁, λ₂, ..., λₙ be distinct eigenvalues of n×n matrix [A] with |λ₁| > |λ₂| > ... > |λₙ|, and v₁, v₂, ..., vₙ be corresponding eigenvectors. Any vector v in the space spanned by eigenvectors can be written as v = c₁v₁ + c₂v₂ + ... + cₙvₙ.

Pre-multiplying by A: Av = λ₁(c₁v₁ + c₂(λ₂/λ₁)v₂ + ... + cₙ(λₙ/λ₁)vₙ)

After r iterations: Aʳv = λ₁ʳ[c₁v₁ + c₂(λ₂/λ₁)ʳv₂ + ... + cₙ(λₙ/λ₁)ʳvₙ]

The eigenvalue λ₁ can be computed as the limit of the ratio of corresponding components of Aʳ⁺¹v and Aʳv: λ₁ = lim(r→∞) (Aʳ⁺¹v)ₚ/(Aʳv)ₚ for any component p.

Finding the Smallest Eigenvalue

For finding the eigenvalue of least magnitude, apply the Power Method to [A⁻¹]. Since A⁻¹ = (1/λ)(X), the inverse matrix has eigenvalues that are reciprocals of the eigenvalues of [A].

⭐ Key Takeaways

The Power Method provides an iterative numerical approach for finding the dominant eigenvalue (largest in magnitude) and its associated eigenvector without solving the characteristic polynomial. The method requires choosing an initial vector, repeatedly multiplying by matrix A, and normalizing—the scaling factor at each step approximates the eigenvalue. For convergence, the eigenvalues must be real and distinct, with the dominant eigenvalue being strictly larger in magnitude than all others. If the smallest eigenvalue is needed, the Power Method is applied to the inverse matrix instead.

🧠 Quick Revision Questions

1. What is the core iterative procedure of the Power Method for finding the largest eigenvalue?
2. Why does the ratio of corresponding components of Aʳ⁺¹v and Aʳv converge to the dominant eigenvalue?
3. In the example with matrix [[2,3,2],[4,3,5],[3,2,9]], what was the largest eigenvalue after five iterations?
4. How would you modify the Power Method to find the smallest eigenvalue of a matrix?
5. What condition must the eigenvalues satisfy for the Power Method to converge properly?

📘 Lecture 17 — Jacobi’s Method

📖 Overview: Jacobi’s method is an iterative technique for computing all eigenvalues and eigenvectors of a real symmetric matrix. It works by systematically annihilating off-diagonal elements through a series of orthogonal (rotation) transformations until the matrix becomes diagonal, from which the eigenvalues can be directly read. This method is highly recommended for its numerical stability and ability to produce the complete eigen-system.

🗂️ Topics Covered

The lecture introduces Jacobi’s method for eigenvalue computation of real symmetric matrices, starting with the definition of orthogonal matrices and the underlying diagonalization principle. It details the construction of rotation matrices to zero out the largest off-diagonal element, including the formula for the rotation angle. Three complete worked examples demonstrate the iterative process, including how to compute the eigenvectors from the product of the rotation matrices.

📝 Lecture Summary

Jacobi’s Method

An n×n matrix [A] is said to be orthogonal if [A]ᵀ[A] = [I], meaning [A]ᵀ = [A]⁻¹. For computing all eigenvalues and eigenvectors of a real symmetric matrix, Jacobi’s method is highly recommended. It is based on the property that if [A] is an n×n real symmetric matrix, its eigenvalues are real, and there exists an orthogonal matrix [S] such that the diagonal matrix D = [S⁻¹][A][S]. This diagonalization is carried out by applying a series of orthogonal transformations S₁, S₂, ..., Sₙ.

Let A be an n×n real symmetric matrix. Suppose aᵢⱼ is the numerically largest element among the off-diagonal elements of A. We construct an orthogonal matrix S₁ defined with:

• sᵢᵢ = cosθ, sᵢⱼ = -sinθ
• sⱼᵢ = sinθ, sⱼⱼ = cosθ
• All remaining off-diagonal elements are zero, and remaining diagonal elements are unity.

After the transformation D₁ = S₁⁻¹ A S₁ = S₁ᵀ A S₁, the elements at positions (i,j) and (j,i) are annihilated. The new diagonal elements dᵢᵢ and dⱼⱼ and the off-diagonal element dᵢⱼ are computed from:

• dᵢⱼ = 0 only if aᵢⱼ cos 2θ + (aⱼⱼ - aᵢᵢ) sin 2θ / 2 = 0

🔑 Definition — Orthogonal matrix: A matrix [A] is orthogonal if [A]ᵀ[A] = [I], or equivalently [A]ᵀ = [A]⁻¹. 📐 Formula: tan 2θ = 2aᵢⱼ / (aᵢᵢ - aⱼⱼ) → This determines the rotation angle that zeroes out the off-diagonal pair. 📌 Example: For matrix A = [[1, 2], [3, 2], [2, 1]], the largest off-diagonal element is a₁₃ = a₃₁ = 2. Then tan 2θ = 2(2)/(1-1) = 4/0 = ∞, giving θ = π/4. The orthogonal matrix S₁ = [[1/√2, 0, -1/√2], [0, 1, 0], [1/√2, 0, 1/√2]] annihilates d₁₃ and d₃₁.

However, though it creates a new pair of zeros, it also introduces non-zero contributions at formerly zero positions.

Example 1 — Full Jacobi Iteration

Find all eigenvalues and corresponding eigenvectors of A = [[1, 2], [3, 2], [2, 1]] by Jacobi’s method.

Step 1: Largest off-diagonal element is a₁₃ = a₃₁ = 2. With θ = π/4, S₁ = [[1/√2, 0, -1/√2], [0, 1, 0], [1/√2, 0, 1/√2]]. Then D₁ = S₁⁻¹ A S₁ = [[3, 2, 0], [2, 3, 0], [0, 0, -1]]. Sum of diagonal elements of D₁ equals sum of diagonal elements of A = 3+2+1=6.

Step 2: Largest off-diagonal element of D₁ is d₁₂ = d₂₁ = 2. Again tan 2θ = 2(2)/(3-3) = ∞, giving θ = π/4. S₂ = [[1/√2, -1/√2, 0], [1/√2, 1/√2, 0], [0, 0, 1]]. D₂ = S₂⁻¹ D₁ S₂ = [[5, 0, 0], [0, 1, 0], [0, 0, -1]]. This is diagonal, so eigenvalues are 5, 1, -1.

The eigenvectors are columns of S = S₁S₂ = [[1/√2, -1/√2, -1/√2], [1/√2, 1/√2, 0], [1/√2, -1/√2, 1/√2]].

💡 Why this matters: The sum of diagonal elements (trace) remains invariant under each rotation, providing a check for numerical accuracy.

Example 2 — Multiple Rotations Required

Find all eigenvalues of A = [[2, -1, 0], [-1, 2, -1], [0, -1, 2]] by Jacobi’s method.

Step 1: Choose a₁₂ = -1 as largest element. tan 2θ = 2(-1)/(2-2) = ∞, giving θ = π/4. With cosθ = sinθ = 1/√2, S₁ = [[1/√2, -1/√2, 0], [1/√2, 1/√2, 0], [0, 0, 1]]. D₁ = [[1, 0, -1/√2], [0, 3, -1/√2], [-1/√2, -1/√2, 2]].

Step 2: Largest element is d₁₃ = -1/√2. tan 2θ = 2(-1/√2)/(1-2) = -√2/(-1) = √2, giving θ = 27°22′41″. S₂ = [[0.888, 0, -0.459], [0, 1, 0], [0.459, 0, 0.888]]. D₂ = [[0.634, -0.325, 0], [-0.325, 3, -0.628], [0, -0.628, 2.365]].

Step 3: Largest off-diagonal is d₂₃ = -0.628. tan 2θ = 2(-0.628)/(3-2.365) = -1.256/0.635 = -1.978, giving θ = -31°35′24″. S₃ = [[1, 0, 0], [0, 0.852, 0.524], [0, -0.524, 0.852]]. D₃ = [[0.634, -0.277, 0], [0.277, 3.386, 0], [0, 0, 1.979]].

Further rotations would fully diagonalize, but approximate eigenvalues are 0.634, 3.386, and 1.979.

Example 3 — Real Application

Find eigenvalues and eigenvectors of A = [[1, 1/2, 1/3], [1/2, 1/3, 1/4], [1/3, 1/4, 1/5]] by Jacobi’s method.

Step 1: Largest off-diagonal is a₁₂ = a₂₁ = 1/2. tan 2θ = 2(1/2)/(1 - 1/3) = 1/(2/3) = 3/2 = 1.5. θ = tan⁻¹(1.5)/2 = 28.155°. S₁ = [[0.882, -0.472, 0], [0.472, 0.882, 0], [0, 0, 1]]. D₁ = [[1.268, 0.000, 0.412], [0.000, 0.066, 0.063], [0.412, 0.063, 0.200]]. Sum of diagonals = 1.53, matching original trace.

Step 2: Largest off-diagonal is d₁₃ = d₃₁ = 0.412. tan 2θ = 2(0.412)/(1.268 - 0.200) = 0.824/1.068 = 0.772. θ = tan⁻¹(0.772)/2 = 18.834°. S₂ = [[0.946, 0, -0.323], [0, 1, 0], [0.323, 0, 0.946]]. D₂ = [[1.408, 0.020, -0.001], [0.020, 0.066, 0.060], [-0.001, 0.060, 0.059]]. Sum of diagonals = 1.53.

Eigenvalues are approximately 1.408, 0.066, and 0.059. The eigenvectors are columns of S = S₁S₂ = [[0.8343, -0.472, -0.2848], [0.446, 0.88, -0.1524], [0.323, 0, 0.946]].

⭐ Key Takeaways

Jacobi’s method transforms a real symmetric matrix into diagonal form through successive orthogonal (rotation) transformations, each zeroing the largest off-diagonal element. The rotation angle θ is determined by tan 2θ = 2aᵢⱼ/(aᵢᵢ - aⱼⱼ), which ensures the chosen pair gets annihilated while preserving the matrix's trace. Although each rotation creates new non-zero off-diagonal elements, continued iterations drive all off-diagonal elements toward zero. The eigenvalues appear as diagonal entries of the final matrix, and the eigenvectors are columns of the product of all rotation matrices.

🧠 Quick Revision Questions

1. What condition must the matrix satisfy for Jacobi's method to be applicable?
2. How is the rotation angle θ computed in Jacobi's method?
3. Why does the sum of diagonal elements remain unchanged after each rotation?
4. How are the eigenvectors obtained at the end of the Jacobi iteration?
5. What happens to the off-diagonal elements that were already zeroed when a new rotation annihilates a different pair?

📘 Lecture 18 — Finite Differences

📖 Overview: This lecture introduces the fundamental concept of finite differences, which are crucial for numerical interpolation and extrapolation of tabulated functions. It defines and explores the relationships between forward, backward, and central difference operators, as well as the shift, average, and differential operators, establishing the symbolic foundations necessary for deriving interpolation formulas.

🗂️ Topics Covered

This lecture begins with the definition of interpolation and extrapolation for tabulated functions with equally spaced abscissas. It then introduces the forward difference operator (Δ), shows how to construct a forward difference table, and derives expressions for function values in terms of leading differences. The backward difference operator (∇) and central difference operator (δ) are similarly defined, with their own difference tables. Finally, the shift operator (E), average operator (μ), and differential operator (D) are introduced, along with crucial symbolic relationships linking all these operators (e.g., Δ = E - 1, E = e^(hD)).

📝 Lecture Summary

Introduction

Finite differences are essential for numerical techniques involving tabulated values of functions. For a function (y = f(x)) with given values ((x_k, y_k)), the process of estimating (y) for an intermediate (x) is called interpolation, while estimating (y) for an (x) outside the table is called extrapolation. The concept of finite differences and the symbolic relations between different difference operators are fundamental for establishing various interpolation formulas.

Finite Difference Operators

Forward Differences

For a table of values ((x_k, y_k)) with equally spaced abscissas, the forward difference operator (\Delta) is defined as (\Delta y_i = y_{i+1} - y_i), for (i = 0,1,...,(n-1)). These are called first differences. The differences of first differences are second differences, defined as (\Delta^2 y_i = \Delta y_{i+1} - \Delta y_i). In general, the r-th difference is (\Delta^r y_i = \Delta^{r-1} y_{i+1} - \Delta^{r-1} y_i). These are systematically displayed in a forward difference table (or diagonal difference table), where each difference is placed midway between the elements of the previous column, and the subscript remains constant along each diagonal. The first term (y_0) is the leading term, and the differences (\Delta y_0, \Delta^2 y_0, \Delta^3 y_0,...) are leading differences.

🔑 Definition — Forward Difference Operator (Δ): (\Delta y_i = y_{i+1} - y_i), where (y_i = f(x_i)) for equally spaced (x_i).

📐 Formula: (\Delta^n y_0 = y_n - \ ^nC_1 y_{n-1} + \ ^nC_2 y_{n-2} - \ ^nC_3 y_{n-3} + \cdots + (-1)^n y_0) → The coefficients in the expansion of (\Delta^n y_0) are binomial coefficients with alternating signs. 📌 Example: Express (\Delta^2 y_0) and (\Delta^3 y_0) in terms of the values of the function (y).

• (\Delta^2 y_0 = \Delta y_1 - \Delta y_0 = (y_2 - y_1) - (y_1 - y_0) = y_2 - 2y_1 + y_0)
• (\Delta^3 y_0 = \Delta^2 y_1 - \Delta^2 y_0 = (y_3 - 2y_2 + y_1) - (y_2 - 2y_1 + y_0) = y_3 - 3y_2 + 3y_1 - y_0)

🔑 Definition — Leading Term ((y_0)): The first value in a forward difference table. 🔑 Definition — Leading Differences ((\Delta y_0, \Delta^2 y_0, ...)): The first differences in each column of a forward difference table.

Backward Differences

For a given table of values with equally spaced abscissas, the first backward differences are expressed using the backward difference operator (\nabla) as (\nabla y_i = y_i - y_{i-1}), for (i = n, (n-1),..., 1). The second backward differences are (\nabla^2 y_i = \nabla y_i - \nabla y_{i-1}), and in general, the k-th backward differences are (\nabla^k y_i = \nabla^{k-1} y_i - \nabla^{k-1} y_{i-1}). These are arranged in a backward difference table, where the subscript remains constant along every backward diagonal. Any value of y can be expressed in terms of (y_n) and its backward differences, leading to symbolic relations like (y_{n-1} = (1-\nabla) y_n).

🔑 Definition — Backward Difference Operator (∇): (\nabla y_i = y_i - y_{i-1}), where (y_i = f(x_i)) for equally spaced (x_i). 📐 Formula: (y_{n-r} = y_n - \ ^nC_1 \nabla y_n + \ ^nC_2 \nabla^2 y_n - \cdots + (-1)^r \nabla^r y_n) → The value of y at a point can be expressed in terms of a base value and its backward differences. 📌 Example: Show that any value of y can be expressed in terms of (y_n) and its backward differences.

• From (\nabla y_i = y_i - y_{i-1}), we get (y_{n-1} = y_n - \nabla y_n).
• From (\nabla^2 y_i = \nabla y_i - \nabla y_{i-1}), we get (\nabla y_{n-1} = \nabla y_n - \nabla^2 y_n).
• Then (y_{n-2} = y_{n-1} - \nabla y_{n-1} = (y_n - \nabla y_n) - (\nabla y_n - \nabla^2 y_n) = y_n - 2\nabla y_n + \nabla^2 y_n). Symbolically, (y_{n-2} = (1 - \nabla)^2 y_n).

Central Differences

In some applications, the central difference notation is more convenient. The central difference operator (\delta) is defined by (\delta y_i = y_{i+(1/2)} - y_{i-(1/2)}). Higher order central differences are defined as (\delta^2 y_i = \delta y_{i+(1/2)} - \delta y_{i-(1/2)}) and (\delta^n y_i = \delta^{n-1} y_{i+(1/2)} - \delta^{n-1} y_{i-(1/2)}). In a central difference table, odd differences have fractional suffixes, and even differences with the same subscript lie horizontally.

🔑 Definition — Central Difference Operator (δ): (\delta y_i = y_{i+(1/2)} - y_{i-(1/2)}), where the subscript is the average of the subscripts of the two y-values.

Shift Operator, E

The shift operator (E) is defined by the property (E f(x) = f(x + h)). Applying it n times gives (E^n f(x) = f(x + nh)). The inverse operator is (E^{-1} f(x) = f(x - h)). This operator is fundamental for deriving symbolic relations between difference operators.

🔑 Definition — Shift Operator (E): (E f(x) = f(x + h)), where (h) is the spacing between consecutive x-values. 📐 Formula: (E = e^{hD}) → The shift operator is equivalent to the exponential of the differential operator times the step size. 📌 Example: Prove that (hD = \log E).

• This stems from the Taylor series expansion: (E y_x = y_{x+h} = y_x + h y'_x + \frac{h^2}{2!} y''_x + \cdots = (1 + hD + \frac{h^2 D^2}{2!} + \cdots) y_x = e^{hD} y_x). Therefore, (E = e^{hD}) and taking logs gives (hD = \log E).

Average Operator, μ

The average operator (\mu) is defined as (\mu f(x) = \frac{1}{2} \left[ f(x + \frac{h}{2}) + f(x - \frac{h}{2}) \right]). It gives the average of the function values at two symmetrically placed points.

🔑 Definition — Average Operator (μ): (\mu f(x) = \frac{1}{2} \left[ f(x + \frac{h}{2}) + f(x - \frac{h}{2}) \right]). 📐 Formula: (\mu = \frac{1}{2}(E^{1/2} + E^{-1/2})) → This shows the relationship between the average operator and the shift operator.

Differential Operator, D

The differential operator (D) is defined as (D f(x) = \frac{d}{dx} f(x) = f'(x)), with higher powers representing higher-order derivatives, e.g., (D^2 f(x) = f''(x)).

Important Results Using {Δ, ∇, δ, E, μ}

By defining the operators on a common variable (y_x), the following fundamental relationships are derived:

• (\Delta = E - 1)
• (\nabla = 1 - E^{-1} = \frac{E-1}{E})
• (\delta = E^{1/2} - E^{-1/2})
• (\mu = \frac{1}{2}(E^{1/2} + E^{-1/2}))
• (E = e^{hD})
• Therefore, (hD = \log E = \log(1 + \Delta) = -\log(1 - \nabla) = \sinh^{-1}(\mu \delta)) 💡 Why this matters: These symbolic relations allow one to express any difference operator in terms of another, which is key to deriving various interpolation formulas and converting between them.

📌 Example: Prove the following results:

• (1) (1 + \mu^2 \delta^2 = \left(1 + \frac{\delta^2}{2}\right)^2)
• (2) (E^{1/2} = \mu + \frac{\delta}{2})
• (3) (\Delta = \frac{\delta^2}{2} + \delta \sqrt{1 + \frac{\delta^2}{4}})
• (4) (\mu \delta = \frac{\Delta}{2} + \frac{\Delta}{2E})
• (5) (\mu \delta = \frac{1}{2}(\Delta + \nabla))
• Solution (1): (\mu \delta = \frac{1}{2}(E^{1/2} + E^{-1/2})(E^{1/2} - E^{-1/2}) = \frac{1}{2}(E - E^{-1})). Then (1 + \mu^2 \delta^2 = 1 + \frac{1}{4}(E^2 - 2 + E^{-2}) = \frac{1}{4}(E + E^{-1})^2). Also, (1 + \frac{\delta^2}{2} = 1 + \frac{1}{2}(E^{1/2} - E^{-1/2})^2 = \frac{1}{2}(E + E^{-1})). Therefore, (\left(1 + \frac{\delta^2}{2}\right)^2 = \frac{1}{4}(E + E^{-1})^2 = 1 + \mu^2 \delta^2).
• Solution (2): (\mu + \frac{\delta}{2} = \frac{1}{2}(E^{1/2} + E^{-1/2}) + \frac{1}{2}(E^{1/2} - E^{-1/2}) = E^{1/2}).
• Solution (3): From (\Delta = E - 1) and the proof in (2), we have (E^{1/2} = \mu + \frac{\delta}{2}). Noting that (\mu = \sqrt{1 + \frac{\delta^2}{4}}), we get (\Delta = E - 1 = (E^{1/2})^2 - 1 = (\mu + \frac{\delta}{2})^2 - 1 = \mu^2 + \mu \delta + \frac{\delta^2}{4} - 1 = (1 + \frac{\delta^2}{4}) + \mu \delta + \frac{\delta^2}{4} - 1 = \frac{\delta^2}{2} + \mu \delta). Since (\mu = \sqrt{1 + \frac{\delta^2}{4}}), we have (\mu \delta = \delta \sqrt{1 + \frac{\delta^2}{4}}), so (\Delta = \frac{\delta^2}{2} + \delta \sqrt{1 + \frac{\delta^2}{4}}).
• Solution (4): (\mu \delta = \frac{1}{2}(E - E^{-1}) = \frac{1}{2}(1 + \Delta - E^{-1}) = \frac{\Delta}{2} + \frac{1}{2}(1 - E^{-1}) = \frac{\Delta}{2} + \frac{1}{2}\frac{E-1}{E} = \frac{\Delta}{2} + \frac{\Delta}{2E}).
• Solution (5): (\mu \delta = \frac{1}{2}(E - E^{-1}) = \frac{1}{2}(1 + \Delta - 1 + \nabla) = \frac{1}{2}(\Delta + \nabla)).

⭐ Key Takeaways

A student must understand the definitions and graphical representation of forward (Δ), backward (∇), and central (δ) finite difference operators for equally spaced data. The ability to derive symbolic relationships between these operators, particularly using the shift operator (E) as a bridge, is absolutely critical. The key formulas to memorize are Δ = E-1, ∇ = 1-E⁻¹, δ = E^(1/2) - E^(-1/2), μ = (E^(1/2)+E^(-1/2))/2, and E = e^(hD). From these, one can derive that functions of one operator (like log(1+Δ), sinh⁻¹(μδ)) can be expressed in terms of others, which is the foundation for building different interpolation formulas. Finally, the skill of expressing a function value in terms of its leading differences (e.g., (y_n = (1+\Delta)^n y_0) and (y_{n-r} = (1-\nabla)^r y_n)) is a core technique.

🧠 Quick Revision Questions

1. Define the forward difference operator Δ and the backward difference operator ∇ for a function y = f(x) with equally spaced abscissas.
2. Show that (\Delta^3 y_0 = y_3 - 3y_2 + 3y_1 - y_0).
3. Using the shift operator E, derive the symbolic relationship (\Delta) = E - 1.
4. Prove that the differential operator D and the shift operator E are related by the formula E = e^{hD}.
5. Derive the expression for the average operator μ in terms of the shift operator E.

📘 Lecture 19 — Interpolation Newton’s Forward Difference Formula

📖 Overview: This lecture introduces Newton’s Forward Difference Formula, a powerful method for interpolating function values at equally spaced points. The formula expresses the value at any point using the initial function value and its forward differences, making it especially useful for estimating values near the beginning of a data table.

🗂️ Topics Covered

The lecture begins with the derivation of Newton’s Forward Difference Formula using the shift operator E and forward difference operator Δ. It then presents the complete formula with binomial coefficients, explains its alternative form using y notation, and discusses its application for interpolation and extrapolation near the start of tabular data. Several worked examples demonstrate the formula’s use for evaluating specific function values, finding interpolating polynomials, and estimating missing table entries.

📝 Lecture Summary

Introduction and Formula Derivation

Let y = f(x) be a function that takes values at equally spaced points x₀, x₀ + h, x₀ + 2h, ... with spacing h. To evaluate f(x₀ + ph) for any real number p, we define the shift operator E such that Eᵖ f(x) = f(x + ph). Using the relationship E = 1 + Δ, we expand (1 + Δ)ᵖ as a binomial series to obtain Newton’s Forward Difference Formula.

The formula expresses f(x₀ + ph) in terms of f(x₀) and its leading forward differences:

f(x₀ + ph) = f(x₀) + pΔf(x₀) + [p(p-1)/2!] Δ²f(x₀) + [p(p-1)(p-2)/3!] Δ³f(x₀) + ... + [p(p-1)...(p-n+1)/n!] Δⁿf(x₀) + Error

An alternative expression uses y notation: yₓ = y₀ + pΔy₀ + [p(p-1)/2!] Δ²y₀ + [p(p-1)(p-2)/3!] Δ³y₀ + ... + [p(p-1)...(p-n+1)/n!] Δⁿy₀ + Error

Here, p = (x - x₀)/h. If we retain (r+1) terms, we obtain a polynomial of degree r agreeing with y at points x₀, x₁, ..., xᵣ.

🔑 Definition — Newton’s Forward Difference Formula: A formula for interpolating function values at equally spaced points using forward differences, primarily used for interpolation near the beginning of a set of tabular values and for extrapolation a short distance backward from y₀.

📐 Formula: f(x₀ + ph) = f(x₀) + pΔf(x₀) + [p(p-1)/2!] Δ²f(x₀) + [p(p-1)(p-2)/3!] Δ³f(x₀) + ... → This gives the function value at any point x = x₀ + ph by using the initial value and successive forward differences.

💡 Why this matters: This formula is fundamental for numerical interpolation; it allows us to estimate function values between known data points when the data is evenly spaced.

Example 1: Evaluating f(15)

Given the table of values: f(10)=46, f(20)=66, f(30)=81, f(40)=93, f(50)=101.

First, construct the forward difference table:

• Δy₀ = 20, Δ²y₀ = -5, Δ³y₀ = 2, Δ⁴y₀ = -3

Here, x₀ = 10, y₀ = 46, h = 10, and p = (15 - 10)/10 = 0.5

Using Newton’s Forward Difference Formula: f(15) = 46 + 0.5(20) + 0.5(0.5-1)/2 + 0.5(0.5-1)(0.5-2)/6 + 0.5(0.5-1)(0.5-2)(0.5-3)/24

f(15) = 46 + 10 + 0.625 + 0.125 + 0.1172 = 56.8672

So f(15) = 56.8672 correct to four decimal places.

📌 Example: Using the formula with p=0.5 gave f(15)=56.8672, where each term contributed: 46 + 10 + (-0.625) correction + 0.125 + 0.1172.

Example 2: Finding Newton’s Forward Difference Interpolating Polynomial

Given data: x=0.1,0.2,0.3,0.4; y=1.40,1.56,1.76,2.00

Construct the forward difference table:

• Δy₀ = 0.16, Δ²y₀ = 0.04 (third and fourth leading differences are zero)

Since Δ³y₀ = 0, we have: y = y₀ + pΔy₀ + [p(p-1)/2] Δ²y₀

Here x₀ = 0.1, h = 0.1, p = (x - 0.1)/0.1 = 10x - 1

Substituting: y = 1.40 + (10x - 1)(0.16) + (10x - 1)(10x - 2)/2

This is the required Newton’s interpolating polynomial.

📌 Example: When differences beyond the second are zero, the formula reduces to a quadratic polynomial in terms of x.

Example 3: Estimating the Missing Figure

Given table: x=0,1,2,4; y=2,5,7,32 (missing at x=3)

Since we have four entries, the function can be represented by a polynomial of degree three, so Δ⁴f(x) = 0 for all x.

Using Δ⁴f(x₀) = 0, which means (E - 1)⁴ f(x₀) = 0, expanding: E⁴ - 4E³ + 6E² - 4E + 1 = 0

This gives: f(x₄) - 4f(x₃) + 6f(x₂) - 4f(x₁) + f(x₀) = 0

Using the given values: 32 - 4f(x₃) + 6(7) - 4(5) + 2 = 0 → 32 - 4f(x₃) + 42 - 20 + 2 = 0 → -4f(x₃) = -56 → f(x₃) = 14

The missing value at x=3 is 14.

📌 Example: By using the property that fourth differences are zero for a cubic polynomial, we derived a linear equation to solve for the missing value.

Example 4: Finding f(0.36)

Given table: x=0.2,0.3,0.4,0.5,0.6; y=0.2304,0.2788,0.3222,0.3617,0.3979

Forward difference table:

• Δy₀ = 0.0484, Δ²y₀ = -0.005, Δ³y₀ = 0.0011, Δ⁴y₀ = -0.0005

Here x₀ = 0.2, y₀ = 0.2304, h = 0.1, p = (0.36 - 0.2)/0.1 = 1.6

Using formula up to fourth difference: y = 0.2304 + 1.6(0.0484) + 1.6(1.6-1)/2 + 1.6(1.6-1)(1.6-2)/6 + 1.6(1.6-1)(1.6-2)(1.6-3)/24

y = 0.2304 + 0.07744 - 0.0024 - 0.00007 - 0.00001 = 0.3053

So f(0.36) = 0.3053

📌 Example: With p=1.6, each term was calculated: 0.2304 + 0.077441 - 0.0024 - 0.00007 - 0.00001 = 0.3053.

Example 5: Finding a Cubic Polynomial

Given data: x=0,1,2,3,4,5; y=-3,3,11,27,57,107

Construct the difference table: Δy₀ = 6, Δ²y₀ = 2, Δ³y₀ = 6 (fourth and higher differences are zero)

Here h = 1, p = x - 0 = x, x₀ = 0

Using Newton’s Forward Difference Formula up to third difference: f(x) = -3 + 6x + x(x-1)/2 + x(x-1)(x-2)/6

f(x) = -3 + 6x + x(x-1) + x(x-1)(x-2)

f(x) = -3 + 6x + x² - x + x³ - 3x² + 2x

f(x) = x³ - 2x² + 7x - 3

This is the required cubic polynomial.

📌 Example: Since third differences were constant (6) and fourth differences zero, the polynomial was cubic, and substituting p=x gave the interpolating polynomial.

⭐ Key Takeaways

Newton’s Forward Difference Formula expresses f(x₀ + ph) = f(x₀) + pΔf(x₀) + [p(p-1)/2!]Δ²f(x₀) + [p(p-1)(p-2)/3!]Δ³f(x₀) + ... where p = (x - x₀)/h. This formula is primarily used for interpolation near the beginning of tabular data; it works only for equally spaced points. When higher-order differences become zero, the polynomial degree reduces accordingly. The formula is derived from the binomial expansion of (1+Δ)ᵖ, and it can be used to find missing values by setting appropriate higher differences to zero.

🧠 Quick Revision Questions

1. What is the value of p in Newton’s Forward Difference Formula when x = x₀ + 2.5h?
2. For a cubic polynomial, which forward differences are zero?
3. How is the shift operator E related to the forward difference operator Δ?
4. When interpolating at x = x₀ + 0.3h, what does p equal?
5. If Δ⁴y₀ = 0 for all tabular points, what is the degree of the interpolating polynomial?

📘 Lecture 20 — Newton’s Backward Difference Interpolation Formula

📖 Overview: This lecture introduces Newton’s backward difference interpolation formula, which is used for interpolating the value of a function near the end of a table of values and for extrapolating a short distance forward. It is derived using the backward difference operator and is also known as the Newton-Gregory backward difference interpolation formula. The lecture includes the derivation of the formula and several worked examples demonstrating its application.

🗂️ Topics Covered

The lecture covers the definition and purpose of Newton’s backward difference interpolation formula, its derivation using the shift operator E and the backward difference operator ∇, the binomial expansion that yields the formula, and the expression for the parameter p. It then provides three detailed examples that show how to construct backward difference tables and apply the formula to estimate function values at specific points, such as f(7.5), sales for the year 1979, and f(1.45).

📝 Lecture Summary

Newton’s Backward Difference Interpolation Formula

The Newton’s backward difference interpolation formula is used to find the value of a function (y = f(x)) near the end of a table of values or to extrapolate a value forward. Let the function take values (f(x_n), f(x_n - h), f(x_n - 2h), \ldots, f(x_0)) at equally spaced points (x_n, x_n - h, x_n - 2h, \ldots, x_0). To evaluate the function at ((x_n + ph)), where p is any real number, the shift operator (E) is used. The derivation begins with the relation (f(x_n + ph) = E^p f(x_n) = (E^{-1})^{-p} f(x_n) = (1 - \nabla)^{-p} f(x_n)). Applying the binomial expansion yields the formula.

🔑 Definition — 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)\cdots(p+n-1)}{n!} \nabla^n f(x_n) + \text{Error}). This formula is also known as the Newton-Gregory backward difference interpolation formula. If (r+1) terms are retained, we obtain a polynomial of degree (r) agreeing with (f(x)) at (x_n, x_{n-1}, \ldots, x_{n-r}).

📐 Formula: (y_x = y_n + p\nabla y_n + \frac{p(p+1)}{2!} \nabla^2 y_n + \frac{p(p+1)(p+2)}{3!} \nabla^3 y_n + \cdots + \frac{p(p+1)\cdots(p+n-1)}{n!} \nabla^n y_n + \text{Error}), where (p = \frac{x - x_n}{h}).

Example 1

Problem: For the following table of values, estimate (f(7.5)).

Solution: The value to be interpolated is at the end of the table, so Newton’s backward interpolation formula is appropriate. First, construct the backward difference table.

Since (3^{\text{rd}}) and higher order differences are zero or constant, the required formula is (y_x = y_n + p\nabla y_n + \frac{p(p+1)}{2!} \nabla^2 y_n + \frac{p(p+1)(p+2)}{3!} \nabla^3 y_n). Here, (x_n = 8, y_n = 512, h = 1). So, (p = \frac{x - x_n}{h} = \frac{7.5 - 8.0}{1} = -0.5). From the table, (\nabla y_n = 169, \nabla^2 y_n = 42, \nabla^3 y_n = 6).

📌 Example: (y_{7.5} = 512 + (-0.5)(169) + \frac{(-0.5)(0.5)}{2}(42) + \frac{(-0.5)(0.5)(1.5)}{6}(6)) (y_{7.5} = 512 - 84.5 - 5.25 - 0.375 = 421.875)

Example 2

Problem: The sales for the last five years is given in the table below. Estimate the sales for the year 1979.

Solution: Let the year 1978 correspond to x=0. The backward difference table is:

Here, (x_n = 5, y_n = 57, h = 1). For the year 1979, (x = 1), so (p = \frac{1 - 5}{1} = -4). Using (y_x = y_n + p\nabla y_n + \frac{p(p+1)}{2!} \nabla^2 y_n + \frac{p(p+1)(p+2)}{3!} \nabla^3 y_n + \frac{p(p+1)(p+2)(p+3)}{4!} \nabla^4 y_n), with (\nabla y_n = 5, \nabla^2 y_n = 1, \nabla^3 y_n = 2, \nabla^4 y_n = 5), we get: (y_{1979} = 57 + (-4)(5) + \frac{(-4)(-3)}{2}(1) + \frac{(-4)(-3)(-2)}{6}(2) + \frac{(-4)(-3)(-2)(-1)}{24}(5)) (y_{1979} = 57 - 20 + 6 - 8 + 5 = 40). Therefore, sales for 1979 are estimated to be 40.

💡 Why this matters: This example shows how the formula can be used for extrapolation (estimating a value before the first data point in the table), not just interpolation near the end.

Example 3

Problem: Consider the following table of values. Use Newton’s Backward Difference Formula to estimate the value of (f(1.45)).

Solution: The backward difference table is:

Here, (x_n = 1.5, y_n = 4.5, h = 0.1). For (x = 1.45), (p = \frac{1.45 - 1.5}{0.1} = -0.5). Values from the table: (\nabla y_n = 1, \nabla^2 y_n = 0.2, \nabla^3 y_n = -0.2, \nabla^4 y_n = -0.4, \nabla^5 y_n = -0.5). Using the full formula up to the 5th difference: (y_{1.45} = 4.5 + (-0.5)(1) + \frac{(-0.5)(-0.5+1)}{2!}(0.2) + \frac{(-0.5)(-0.5+1)(-0.5+2)}{3!}(-0.2) + \frac{(-0.5)(-0.5+1)(-0.5+2)(-0.5+3)}{4!}(-0.4) + \frac{(-0.5)(-0.5+1)(-0.5+2)(-0.5+3)(-0.5+4)}{5!}(-0.5)) (y_{1.45} = 4.5 + (-0.5)(1) + \frac{(-0.5)(0.5)}{2}(0.2) + \frac{(-0.5)(0.5)(1.5)}{6}(-0.2) + \frac{(-0.5)(0.5)(1.5)(2.5)}{24}(-0.4) + \frac{(-0.5)(0.5)(1.5)(2.5)(3.5)}{120}(-0.5)) (y_{1.45} = 4.5 - 0.5 - 0.025 + 0.0125 + 0.015625 + 0.068359 = 4.07148)

⭐ Key Takeaways

Newton’s backward difference interpolation formula is the essential tool for estimating function values when the point of interest lies near the end of a tabulated dataset. The parameter (p) must be correctly computed as ((x - x_n)/h), where (x_n) is the last tabulated point. The accuracy of the interpolation depends on the order of differences included, and while the formula uses backward differences, its application mirrors that of its forward counterpart. It is also effective for extrapolating a short distance backward from the table’s last point, as demonstrated in a sales estimation example.

🧠 Quick Revision Questions

1. When should Newton’s backward difference interpolation formula be used instead of the forward difference formula?
2. What is the formula for the parameter (p) in Newton’s backward interpolation formula, and what does each term represent?
3. In Example 1, how does knowing that 3rd order differences are constant simplify the formula used?
4. If you have data points at x = 0, 1, 2, 3 and want to estimate f(2.5), what is the value of (p) assuming you use the backward formula?
5. What is the primary distinction between interpolation and extrapolation in the context of Newton’s backward formula?

📘 Lecture 21 — Lagrange’s Interpolation Formula

📖 Overview: This lecture presents Lagrange’s interpolation formula, which is used to construct an interpolating polynomial when the values of the independent variable (x) are not equally spaced. Unlike Newton’s forward/backward formulas, Lagrange’s method works for any spacing. The lecture derives the formula step-by-step, introduces the Lagrange basis polynomials (L_k(x)), and provides worked examples.

🗂️ Topics Covered

The lecture covers the motivation for Lagrange’s interpolation when data points are not equally spaced, the derivation of the polynomial form using coefficients (a_k), the compact form using Lagrange basis polynomials (L_k(x)) and the Kronecker delta (\delta_{ij}), the alternative product notation using (\Pi(x)) and its derivative, and finally two fully worked numerical examples showing how to compute the polynomial and evaluate it at a given point.

📝 Lecture Summary

Lagrange’s Interpolation Formula Derivation

Newton’s interpolation formulas developed earlier can be used only when the values of the independent variable (x) are equally spaced and the differences of (y) ultimately become small. If the values of (x) are not given at equidistant intervals, we have the basic formula associated with the name of Lagrange.

Let (y = f(x)) be a function that takes the values (y_0, y_1, \ldots, y_n) corresponding to (x_0, x_1, \ldots, x_n). Since there are ((n+1)) values of (y), we can represent the function (f(x)) by a polynomial of degree (n). We write this polynomial in the form:

[ y = f(x) = a_0 (x - x_1)(x - x_2)\cdots(x - x_n) + a_1 (x - x_0)(x - x_2)\cdots(x - x_n) + \cdots + a_n (x - x_0)(x - x_1)\cdots(x - x_{n-1}) ]

Here, the coefficients (a_k) are chosen so that this equation is satisfied by the ((n+1)) pairs ((x_i, y_i)). Substituting (x = x_0) gives:

[ y_0 = a_0 (x_0 - x_1)(x_0 - x_2)\cdots(x_0 - x_n) ]

Therefore:

[ a_0 = \frac{y_0}{(x_0 - x_1)(x_0 - x_2)\cdots(x_0 - x_n)} ]

Similarly, we obtain:

[ a_1 = \frac{y_1}{(x_1 - x_0)(x_1 - x_2)\cdots(x_1 - x_n)} ]

And in general:

[ a_i = \frac{y_i}{(x_i - x_0)(x_i - x_1)\cdots(x_i - x_{i-1})(x_i - x_{i+1})\cdots(x_i - x_n)} ]

And:

[ a_n = \frac{y_n}{(x_n - x_0)(x_n - x_1)\cdots(x_n - x_{n-1})} ]

Substituting the values of (a_0, a_1, \ldots, a_n) back into the polynomial gives Lagrange’s interpolation formula:

[ y = f(x) = \frac{(x - x_1)(x - x_2)\cdots(x - x_n)}{(x_0 - x_1)(x_0 - x_2)\cdots(x_0 - x_n)} y_0 + \frac{(x - x_0)(x - x_2)\cdots(x - x_n)}{(x_1 - x_0)(x_1 - x_2)\cdots(x_1 - x_n)} y_1 + \cdots + \frac{(x - x_0)(x - x_1)\cdots(x - x_{i-1})(x - x_{i+1})\cdots(x - x_n)}{(x_i - x_0)(x_i - x_1)\cdots(x_i - x_{i-1})(x_i - x_{i+1})\cdots(x_i - x_n)} y_i + \cdots ]

This formula can be used whether the values (x_0, x_1, \ldots, x_n) are equally spaced or not.

🔑 Definition — Lagrange Basis Polynomials: The formula can be written in compact form as:

[ y = f(x) = L_0(x) y_0 + L_1(x) y_1 + \cdots + L_n(x) y_n = \sum_{k=0}^n L_k(x) y_k = \sum_{k=0}^n L_k(x) f(x_k) ]

where

[ L_i(x) = \frac{(x - x_0)(x - x_1)\cdots(x - x_{i-1})(x - x_{i+1})\cdots(x - x_n)}{(x_i - x_0)(x_i - x_1)\cdots(x_i - x_{i-1})(x_i - x_{i+1})\cdots(x_i - x_n)} ]

Properties of Lagrange Basis Polynomials and Alternative Form

We can easily observe that (L_i(x_i) = 1) and (L_i(x_j) = 0) for (i \neq j). This is represented using the Kronecker delta notation:

[ L_i(x_j) = \delta_{ij} = \begin{cases} 1, & \text{if } i = j \ 0, & \text{if } i \neq j \end{cases} ]

💡 Why this matters: This property ensures that the polynomial passes exactly through each data point ((x_i, y_i)).

We can also introduce the product notation:

[ \Pi(x) = \prod_{i=0}^n (x - x_i) = (x - x_0)(x - x_1)\cdots(x - x_n) ]

Its derivative (\Pi'(x)) is a sum of ((n+1)) terms, each missing one factor. We define:

[ P_k(x) = \prod_{i \neq k} (x - x_i) ]

which is the same as (\Pi(x)) except the factor ((x - x_k)) is absent. Then:

[ \Pi'(x) = P_0(x) + P_1(x) + \cdots + P_n(x) ]

When (x = x_k), all terms in the above sum vanish except (P_k(x_k)). Hence:

[ \Pi'(x_k) = P_k(x_k) = (x_k - x_0)\cdots(x_k - x_{k-1})(x_k - x_{k+1})\cdots(x_k - x_n) ]

Therefore:

[ L_k(x) = \frac{P_k(x)}{P_k(x_k)} = \frac{P_k(x)}{\Pi'(x_k)} = \frac{\Pi(x)}{(x - x_k) \Pi'(x_k)} ]

Finally, Lagrange’s interpolation polynomial of degree (n) can be written as:

[ y(x) = f(x) = \sum_{k=0}^n \frac{\Pi(x)}{(x - x_k) \Pi'(x_k)} f(x_k) = \sum_{k=0}^n L_k(x) f(x_k) = \sum_{k=0}^n L_k(x) y_k ]

Example 1: Fitting a Polynomial and Finding y(5)

Find Lagrange’s interpolation polynomial fitting the points (y(1) = -3, y(3) = 0, y(4) = 30, y(6) = 132). Hence find (y(5)).

Solution: Using Lagrange’s interpolation formula:

[ y(x) = f(x) = \frac{(x-3)(x-4)(x-6)}{(1-3)(1-4)(1-6)} (-3) + \frac{(x-1)(x-4)(x-6)}{(3-1)(3-4)(3-6)} (0) + \frac{(x-1)(x-3)(x-6)}{(4-1)(4-3)(4-6)} (30) + \frac{(x-1)(x-3)(x-4)}{(6-1)(6-3)(6-4)} (132) ]

On simplification, we get:

[ y(x) = \frac{1}{10}(-5x^3 + 135x^2 - 460x + 300) = \frac{1}{2}(-x^3 + 27x^2 - 92x + 60) ]

This is the required Lagrange’s interpolation polynomial. Now, (y(5) = 75).

Example 2: Evaluating f(3)

Given the following data, evaluate (f(3)) using Lagrange’s interpolating polynomial.

Solution: Using Lagrange’s formula:

[ f(3) = \frac{(3-2)(3-5)}{(1-2)(1-5)}(1) + \frac{(3-1)(3-5)}{(2-1)(2-5)}(4) + \frac{(3-1)(3-2)}{(5-1)(5-2)}(10) ]

[ f(3) = \frac{(1)(-2)}{(-1)(-4)}(1) + \frac{(2)(-2)}{(1)(-3)}(4) + \frac{(2)(1)}{(4)(3)}(10) = \frac{-2}{4}(1) + \frac{-4}{-3}(4) + \frac{2}{12}(10) ]

[ f(3) = -0.5 + \frac{16}{3} + \frac{20}{12} = -0.5 + 5.3333 + 1.6667 = 6.49999 \approx 6.5 ]

Example 3: Finding the Interpolating Polynomial

Find the interpolating polynomial for the data using Lagrange’s formula:

Solution: Using Lagrange’s formula:

[ y(x) = \frac{(x-2)(x+4)}{(1-2)(1+4)}(3) + \frac{(x-1)(x+4)}{(2-1)(2+4)}(-5) + \frac{(x-1)(x-2)}{(-4-1)(-4-2)}(-4) ]

[ y(x) = -\frac{3}{5}(x^2+2x-8) - \frac{5}{6}(x^2+3x-4) - \frac{4}{30}(x^2-3x+2) ]

[ y(x) = -\frac{3}{5}x^2 - \frac{6}{5}x + \frac{24}{5} - \frac{5}{6}x^2 - \frac{15}{6}x + \frac{20}{6} - \frac{4}{30}x^2 + \frac{12}{30}x - \frac{8}{30} ]

Combining coefficients:

[ y(x) = \left(-\frac{3}{5} - \frac{5}{6} - \frac{4}{30}\right)x^2 + \left(-\frac{6}{5} - \frac{15}{6} + \frac{12}{30}\right)x + \left(\frac{24}{5} + \frac{20}{6} - \frac{8}{30}\right) ]

[ y(x) = -\frac{47}{30}x^2 - \frac{33}{10}x + \frac{118}{15} ]

Which is the required polynomial.

⭐ Key Takeaways

Lagrange’s interpolation formula is a fundamental method for constructing an interpolating polynomial that works for any set of data points, regardless of whether the (x)-values are equally spaced. The key to the formula is the set of Lagrange basis polynomials (L_k(x)), which are designed to be exactly 1 at (x = x_k) and 0 at all other data points, ensuring the polynomial passes through all given points. The formula is completely general for ((n+1)) data points producing a polynomial of degree (n). For the exam, remember the compact sum form (y(x) = \sum L_k(x) y_k), the Kronecker delta property ((L_i(x_j) = \delta_{ij})), and be able to apply the formula step-by-step to small datasets (3-4 points) as shown in the examples.

🧠 Quick Revision Questions

1. What is the main advantage of Lagrange’s interpolation over Newton’s forward/backward formulas?
2. Write down the general expression for the Lagrange basis polynomial (L_i(x)).
3. What is the value of (L_i(x_j)) when (i = j) and when (i \neq j)? What is this property called?
4. In the example with points ((1,-3), (3,0), (4,30), (6,132)), what is the value of (y(5))?
5. If you have 5 data points, what is the maximum degree of the Lagrange interpolating polynomial that fits them?

📘 Lecture 22 — Divided Differences

📖 Overview: This lecture introduces the concept of divided differences, a fundamental tool for constructing interpolation polynomials when data points are not equally spaced. It defines divided differences recursively, presents the standard symmetric form, and develops Newton's divided difference interpolation formula, showing its equivalence to Lagrange's form while requiring fewer arithmetic operations. The lecture also covers the error term associated with this interpolation method.

🗂️ Topics Covered

The lecture begins by defining zero-th, first, and higher-order divided differences recursively. It then presents the standard symmetric form of divided differences and provides an example of constructing a Newton's divided difference table. Newton's divided difference interpolation formula is derived by choosing coefficients that satisfy the data points, leading to a lower triangular system. The lecture shows that Newton's formula can be expressed in terms of forward, backward, and central differences for equally spaced data. It demonstrates through examples that Newton's and Lagrange's interpolation formulas yield the same polynomial, with Newton's formula being computationally more efficient. Finally, the error term for interpolation is derived using the (n+1)th derivative of the function.

📝 Lecture Summary

DIVIDED DIFFERENCES

Assume the function y = f(x) is known for several values of x, (xᵢ, yᵢ), for i = 0, 1, ..., n. The divided differences of orders 0, 1, 2, ..., n are defined recursively. The zero-th order divided difference is y[x₀] = y(x₀) = y₀. The first order divided difference is defined as y[x₀, x₁] = (y₁ - y₀) / (x₁ - x₀). Similarly, higher order divided differences are defined in terms of lower order divided differences by the relation y[x₀, x₁, ..., xₙ] = (y[x₁, x₂, ..., xₙ] - y[x₀, x₁, ..., xₙ₋₁]) / (xₙ - x₀).

We can easily verify that the divided difference is a symmetric function of its arguments. That is, y[x₀, x₁] = y[x₁, x₀] = y₀/(x₀ - x₁) + y₁/(x₁ - x₀). The general symmetric result is y[x₀, ..., xₖ] = Σᵢ₌₀ᵏ [yᵢ / Πⱼ₌₀,ⱼ≠ᵢᵏ (xᵢ - xⱼ)].

📌 Example: Construct the Newton's divided difference table for values of x = 1, 2, 3, 4, 5, 6 and f(x) = -3, 0, 15, 48, 105, 192.

Solution: The table is:

🔑 Definition — First order divided difference: y[x₀, x₁] = (y₁ - y₀) / (x₁ - x₀) 📐 Formula: Divided difference recursion: y[x₀, x₁, ..., xₙ] = (y[x₁, x₂, ..., xₙ] - y[x₀, x₁, ..., xₙ₋₁]) / (xₙ - x₀) → The nth order divided difference is the difference of two (n-1)th order differences divided by the interval between the first and last x-values. 📐 Formula: Symmetric form: y[x₀, ..., xₖ] = Σᵢ₌₀ᵏ [yᵢ / Πⱼ₌₀,ⱼ≠ᵢᵏ (xᵢ - xⱼ)] → The divided difference can be expressed as a sum of the function values each divided by the product of differences from all other points.

NEWTON'S DIVIDED DIFFERENCE INTERPOLATION FORMULA

Let y = f(x) be a function which takes values y₀, y₁, ..., yₙ corresponding to x = xᵢ, i = 0, 1, ..., n. We choose an interpolating polynomial in the form y = f(x) = a₀ + a₁(x - x₀) + a₂(x - x₀)(x - x₁) + ... + aₙ(x - x₀)(x - x₁)...(x - xₙ₋₁). The coefficients aₖ are chosen so as to satisfy the (n+1) pairs (xᵢ, yᵢ). This gives a lower triangular system. The first equation gives a₀ = y(x₀) = y₀. The second gives a₁ = (y₁ - y₀) / (x₁ - x₀) = y[x₀, x₁]. The third yields a₂ = y[x₀, x₁, x₂]. Thus, Newton's divided difference interpolation formula can be written as y = f(x) = y₀ + (x - x₀) y[x₀, x₁] + (x - x₀)(x - x₁) y[x₀, x₁, x₂] + ... + (x - x₀)(x - x₁)...(x - xₙ₋₁) y[x₀, x₁, ..., xₙ].

Newton's divided differences can also be expressed in terms of forward, backward and central differences. Assuming equally spaced values of abscissa with spacing h, we have y[x₀, x₁] = Δy₀/h, y[x₀, x₁, x₂] = Δ²y₀/(2!h²), and in general y[x₀, x₁, ..., xₙ] = Δⁿy₀/(n!hⁿ). Similarly, in terms of backward differences, y[x₀, x₁] = ∇y₁/h, y[x₀, x₁, x₂] = ∇²y₂/(2!h²), and in general y[x₀, x₁, ..., xₙ] = ∇ⁿyₙ/(n!hⁿ). In terms of central differences, y[x₀, x₁] = δy₁/₂/h, y[x₀, x₁, x₂] = δ²y₁/(2!h²), and for even and odd orders: y[x₀, x₁, ..., x₂ₘ] = δ²ᵐyₘ/((2m)!h²ᵐ) and y[x₀, x₁, ..., x₂ₘ₊₁] = δ²ᵐ⁺¹yₘ₊₁/₂/((2m+1)!h²ᵐ⁺¹).

📌 Example: Find the interpolating polynomial by (i) Lagrange's formula and (ii) Newton's divided difference formula for the data: x = 0, 1, 2, 4 and y = 1, 1, 2, 5. Hence show they represent the same polynomial.

Solution: The divided difference table is:

(i) Lagrange's formula gives: f(x) = x³/12 - 3x²/4 + 2x/3 + 1 (ii) Newton's formula gives: f(x) = 1 + (x-0)(0) + (x-0)(x-1)(1/2) + (x-0)(x-1)(x-2)(-1/12) = x³/12 - 3x²/4 + 2x/3 + 1 We observe that the interpolating polynomial by both formulae is one and the same.

💡 Why this matters: Newton's formula involves fewer arithmetic operations than Lagrange's.

📌 Example: Using Newton's divided difference formula, find the quadratic equation for the data: x = 0, 1, 4 and y = 2, 1, 4. Hence find y(2).

Solution: The table is:

Newton's formula gives: y = 2 + (x-0)(-1) + (x-0)(x-1)(1/2) = (x² - 3x + 4)/2. Hence y(2) = 1.

📌 Example: Find the equation of a cubic curve through points (4, -43), (7, 83), (9, 327), (12, 1053) using Newton's divided difference formula.

Solution: The table is:

Newton's formula: f(x) = -43 + (x-4)(42) + (x-4)(x-7)(16) + (x-4)(x-7)(x-9)(1) = x³ - 4x² - 7x - 15.

📌 Example: Show that Newton's divided difference and Lagrange's interpolation formulas for three points are identical.

Solution: The Newton formula is y = y₀ + (x-x₀)y[x₀,x₁] + (x-x₀)(x-x₁)y[x₀,x₁,x₂]. Substituting the definitions of divided differences and simplifying yields y = [(x-x₁)(x-x₂)/((x₀-x₁)(x₀-x₂))]y₀ + [(x-x₀)(x-x₂)/((x₁-x₀)(x₁-x₂))]y₁ + [(x-x₀)(x-x₁)/((x₂-x₀)(x₂-x₁))]y₂, which is Lagrange's form.

NEWTON'S DIVIDED DIFFERENCE FORMULA WITH ERROR TERM

Following the basic definition of divided differences, we have the recursive expressions: y(x) = y₀ + (x-x₀)y[x,x₀]; y[x,x₀] = y[x₀,x₁] + (x-x₁)y[x,x₀,x₁]; ... y[x,x₀,...,xₙ₋₁] = y[x₀,x₁,...,xₙ] + (x-xₙ)y[x,x₀,...,xₙ]. Multiplying the second equation by (x-x₀), the third by (x-x₀)(x-x₁), etc., and adding, we obtain y(x) = y₀ + (x-x₀)y[x₀,x₁] + (x-x₀)(x-x₁)y[x₀,x₁,x₂] + ... + (x-x₀)(x-x₁)...(x-xₙ₋₁)y[x₀,x₁,...,xₙ] + ε(x), where ε(x) = (x-x₀)(x-x₁)...(x-xₙ)y[x,x₀,...,xₙ]. The error term vanishes for x = x₀, x₁, ..., xₙ.

ERROR TERM IN INTERPOLATION FORMULAE

If y(x) is approximated by a polynomial Pₙ(x) of degree n, the error is given by ε(x) = y(x) - Pₙ(x) = (x-x₀)(x-x₁)...(x-xₙ)y[x,x₀,...,xₙ]. Alternatively, it can be expressed as ε(x) = ∏(x) y[x,x₀,...,xₙ] = K∏(x). By considering a function F(x) = y(x) - Pₙ(x) - K∏(x) that vanishes at x₀, x₁, ..., xₙ and an arbitrary point x, and applying Rolle's theorem repeatedly, it can be shown that y⁽ⁿ⁺¹⁾(ξ) appears. Since Pₙ(x) is of degree n, its (n+1)th derivative is zero, and ∏⁽ⁿ⁺¹⁾(x) = (n+1)!. Therefore K = y⁽ⁿ⁺¹⁾(ξ) / (n+1)! for some ξ in the interval.

🔑 Definition — Error term: ε(x) = ∏(x) y⁽ⁿ⁺¹⁾(ξ) / (n+1)! 📐 Formula: ε(x) = (x-x₀)(x-x₁)...(x-xₙ) y⁽ⁿ⁺¹⁾(ξ) / (n+1)! → The error in replacing y(x) by an nth degree interpolation polynomial is the product of the differences between x and each data point, multiplied by the (n+1)th derivative of y evaluated at some point ξ, divided by (n+1)!.

⭐ Key Takeaways

Divided differences provide a systematic way to construct interpolation polynomials for unequally spaced data, with the nth order divided difference defined recursively as the difference of two (n-1)th order differences divided by the interval between the extreme x-values. Newton's divided difference interpolation formula uses these differences as coefficients in a polynomial expressed as progressively nested factors (x-xᵢ), making it computationally more efficient than Lagrange's formula since adding a new data point only requires computing one additional term. The divided difference is a symmetric function of its arguments, and for equally spaced data it reduces to forward, backward, or central differences scaled by factorial and step size terms. The error in interpolation is proportional to the product of distances from interpolation points and the (n+1)th derivative, highlighting that higher-degree polynomials can be more accurate for smooth functions but potentially worse near boundaries.

🧠 Quick Revision Questions

1. Define the first-order divided difference y[x₀, x₁] and write its formula in terms of y₀ and y₁.
2. Write Newton's divided difference interpolation formula for n data points (x₀, y₀), (x₁, y₁), ..., (xₙ, yₙ).
3. For equally spaced data with step size h, what is y[x₀, x₁, x₂] in terms of forward differences?
4. State the formula for the error term ε(x) when a function y(x) is approximated by an nth degree interpolation polynomial.
5. Show that for three points (x₀, y₀), (x₁, y₁), (x₂, y₂), Newton's and Lagrange's interpolation formulas are equivalent.