#11 (07/14/2025)

Roots of f(x) = 0

  1. Algebraic equations (Example: x5 + 3 x2 − 2 x = 0)
  2. Bisection method (slow, steady)
  3. Newton-Raphson method (fast, risk)

  1. Algebraic equations
    According to the Galois theorem, it is not possible to solve general algebraic equations of the fifth order and beyond in closed form.
    1. Quadratic equations
      a x2 + b x + c = 0,     x =
      b ±


      b2−4ac

      2a
      .
      (1)
      Example:
      x2 + 200000 x − 3 = 0.
      (2)
      The exact solutions are
      x1 = −200000,     x2 = 0.000015,
      (3)
      #include <stdio.h>
      #include <math.h>
      int main()
      {
      float a, b, c, disc,x1, x2;
      a=1.0;b=200000;c=-3;
      disc=b*b-4*a*c;
      x1=(-b-sqrt(disc))/(2*a);
      x2=(-b+sqrt(disc))/(2*a);
      printf("x1 = %f, x2 = %f\n",x1,x2); 
      return 0;
      } 
      
      What went wrong ??
      From eq.(1), one obtains
      D
      =
      200002 − 4×1.0 ×(−3)
      (4)
      =
      40000000012,
      (5)
      D
      =
      200000.00003
      (6)
      On the other hand,
      b = 200000.0
      (7)
      Subtraction of two numbers with similar magnitude may result in cancellation error.
      Solution 1 Change float to double
      Solution 2 Noting that
      a x2 + b x +c
      =
      a(xx1) (xx2)
      (8)
      =
      a (x2 − (x1+x2) x + x1 x2)
      (9)
      =
      a x2a (x1+x2) x+ a x1 x2,
      (10)
      so

      x1 x2 = c

      a
      ,
      (11)
      hence
      x2 = c

      x1 a
      = −3

      −200000
      = 0.000015.
      (12)
      Solution 3 Rationalization:

      x2
      =
      b + √D

      2a
      (13)
      =
      (−b−√D)(−b+√D)

      2a (−b − √D)
      (14)
      =
      b2D

      2a (−b − √D)
      ,
      (15)
      so
      x2 = −2000002 − 40000000012

      2 (−200000 − 200000)
      = 0.000015.
      (16)
    2. Cubic equations
      x3 + a x2 + b x+ c = 0.
      (17)
      (Try MathCheat.)
    3. Quartic equations
      x4 + a x3 + b x2+ c x + d = 0.
      (18)
      (Try MathCheat.)
  2. Bisection method
    The mean value theorem:
    If f(x) is continuous over x1xx2 and f′(x) exists over x1 < x < x2 and f(x1) f(x2) < 0, there is at least one point x between x1 and x2 such that f(x)=0.
    bisection1.gif
    1. Choose x1 and x2 such that f(x1) f(x2) < 0.
    2. Set x3 ← (x1+x2)/2.
    3. If f(x1) f(x3) < 0, then, set x2x3.
    4. Else set x1x3.
    5. Until |x1x2| < ϵ or f(x3) = 0, repeat b-d.
    bisection2.gif
    /* Compute the square root of 2 */
    #include <stdio.h>
    #include <math.h>
    #define EPS 1.0e-10
    #define N 100
    
    double f(double x)
    {
    return pow(x,2)-2;
    }
    
    /* start of main */
    int main()
    {
    double x1, x2, x3;
    int count;
    
    printf("Enter xleft and xright separated by space =");
    scanf("%lf %lf", &x1, &x2);
    
    /* sanity check */
    if (f(x1)*f(x2)>0) {printf("Invalid x1 and x2 !\n"); return 0;}
    
    /* bisection start */
    for (count=0;count< N; count++)
     {
     x3= (x1+x2)/2.0;
     if (f(x1)*f(x3)<0 ) x2=x3; else x1=x3;
     if ( f(x3)==0.0 || fabs(x1-x2)< EPS ) break;
    }
    
    printf("iteration = %d\n", count);
    printf("x= %lf\n", x1);
    return 0;
    }
    
    Notes:
    1. fabs is a function to return the absolute value of the argument in math.h (list of functions)
    2. Try plotting a graph by gnuplot first.
    Example: By solving ex=2, obtain approximation of ln 2.
  3. Newton-Raphson Method
    newton1.gif
    The equation of a straight line that passes (a,b) with a slope of m is
    yb = m (xa),
    newton2.gif
    so the tangent line that passes (x1, f(x1)) with the slope of f ′(x1) is

    yf(x1) = f ′(x1) (xx1).
    newton3.gif

    0 − f(x1) = f ′(x1) (xx1),

    x = x1 f(x1)

    f ′(x1)
    ,

    x2 = x1 f(x1)

    f ′(x1)
    .

    x3 = x2 f(x2)

    f ′(x2)
    ,
    General forms:

    xn+1 = xn f(xn)

    f ′(xn)
    .
    Example: √2
    By solving f(x) ≡ x2 − 2 = 0 one can obtain a numerical value of √2.

    f(x)
    =
    x2 −2,
    f ′(x)
    =
    2 x.
    (19)
    Start with x1=2.0 (initial guess), then,
    x1
    =
    2
    (20)
    x2
    =
    2− f(2)

    f ′(2)
    = 2 − 2

    4
    = 1.5
    (21)
    x3
    =
    1.5 − f(1.5)

    f ′(1.5)
    = 1.5 − 0.25

    3
    = 1.41667
    (22)
    x4
    =
    1.4167− f(1.4167)

    f ′(1.4167)
    = … = 1.4142.
    (23)
    The convergence is attained after only 4 iterations.
    Algorithm:
    1. Pick an initial guess, x1.
    2. Make sure f ′(x1) ≠ 0.
    3. Repeat
      1. x2x1f(x1)/f ′(x1)
      2. x1x2
    4. Until |x1x2| ≤ ϵ.
    The Newton-Raphson method fails when f ′(xn) is zero.
    
    
    #include <stdio.h>
    #include <math.h>
    #define EPS 1.0e-10
    
    double f(double x)
    {
    	return x*x-2;
    }
    
    double fp(double x)
    {
    	return 2*x;
    }
    
    double newton(double x)
    {
    	return x - f(x)/fp(x);
    }
    
    int main()
    {
    	double x1, x2;
    	int i;
    
    	printf("Enter initial guess  =");
    	scanf("%lf", &x1);
    
    	if (fp(x1)==0.0) {
    		printf("No convergence.\n");
    		return 0;
    	}
    
    	for (i=0;i<100;i++)
    	{
    		x2=newton(x1);
    		if (fabs(x1-x2)< EPS) break;
    		x1=x2;
    	}
    
    	printf("iteration = %d\n", i);
    	printf("x= %lf\n", x1);
    	return 0;
    }
    
    

    Some interesting examples

    The following examples are useful in emergency when you do not have an access to a computer.
    1. Square root of a ( = √a)
      Let f(x) = x2a, the Newton-Raphson method yields

      xn+1
      =
      xn xn2a

      2 xn
      (24)
      =
      1

      2

      xn + a

      xn

      .
      (25)
      Example (√3 = 1.732 …)
      x1=1,     x2 = 1

      2
      (1 + 3

      1
      ) = 2,    x3 = 1

      2
      (2 + 3

      2
      ) = 1.75,    x4 = 1

      2
      (1.75 + 3

      1.75
      ) = 1.732.
    2. Inverse of a ( = 1/a)
      Let f(x) = 1/xa, the Newton-Raphson method yields
      xn+1
      =
      xn
      1

      xn
      a

      −1

      xn2
      (26)
      =
      xn (2 − a xn).
      (27)
      In case you forgot how to do divisions, you can use this formula.
      Example (1/6 = 0.16667…)

      x1 = 0.2,     x2 = 0.2 (2 − 6×0.2)=0.16,    x3 = 0.16 (2 − 6×0.16) = 0.1664,    x4 = 0.1664 (2 − 6×0.1664) = 0.166666.
    3. 1/√a
      Let f(x) = 1/x2a,

      xn+1
      =
      xn 1/xn2a

      − 2/xn3
      (28)
      =
      xn (3 − a xn2)/2.
      (29)

Summary

Solving simultaneous equations


A x = c,
(30)
or







a11
a12
a1n
a21
a22
a2n
:
:
:
:
an1
an2
ann












x1
x2
:
xn






=





c1
c2
:
cn






,
(31)
or









a11 x1 + a12 x2 + a13 x3 + …+ a1nxn = c1
a21 x1 + a22 x2 + a23 x3 + …+ a2nxn = c2
a31 x1 + a32 x2 + a33 x3 + …+ a3nxn = c3
……
an1 x1 + an2 x2 + an3 x3 + …+ annxn = cn
(32)
For many problems in engineering/science, the size of the matrix A can easily go up to n=1,000,00.

Determinant

  1. 2 ×2 matrices



    a11
    a12
    a21
    a22.



    a11 a22a12 a21
    (33)
  2. 3 ×3 matrices




    a11
    a12
    a13
    a21
    a22
    a23
    a31
    a32
    a33




    a11 a22 a33 + a12 a23 a31+ a21 a32 a13a13 a22a31a12 a21 a33a23 a32 a11
    (34)
  3. The determinant of n ×n (n > 3) matrices can be expressed similar to the above. It consists of n! terms each of which is a product of n variables.

Some properties of determinants

  1. | A B | = | A | | B |
  2. | I | = 1 where I is the identity matrix.
  3. | A |=0 ⇒ A is singular.

Determinant/Cramer's rule

Solving
a11 x1 + a12 x2
=
c1
a21 x1 + a22 x2
=
c2
(35)
yields

x1
=
c1 a22c2 a12

a11 a22a12 a21
=



c1
a12
c2
a22







a11
a12
a21
a22



(36)
x2
=
c2 a11c1 a21

a11 a22a12 a21
=



a11
c1
a21
c2







a11
a12
a21
a22



(37)
Similarly, for three simultaneous equations,

a11 x1 + a12 x2 + a13 x3
=
c1
a21 x1 + a22 x2 + a23 x3
=
c2
a31 x1 + a32 x2 + a33 x3
=
c3
(38)
the solution is expressed as

x1
=




c1
a12
a13
c2
a22
a23
c3
a32
a33









a11
a12
a13
a21
a22
a23
a31
a32
a33




,
(39)
x2
=




a11
c1
a13
a21
c2
a23
a31
c3
a33









a11
a12
a13
a21
a22
a23
a31
a32
a33




,
(40)
x3
=




a11
a12
c1
a21
a22
c2
a31
a32
c3









a11
a12
a13
a21
a22
a23
a31
a32
a33




.
(41)
The above formulas are called Cramer's rule and in general can be extended to higher order matrices. However, it's not practical to use Cramer's rule for a set of more than 4 simultaneous equations. The number of multiplications required for Cramer's rule for n simultaneous equations is n! (n−1) (n+1). For n=10, this amounts to 3,628,800.
Approximate computational time for n simultaneous equations with a 100 MFLOPS computer (an old PC) using Cramer's rule is estimated as shown in Table .1
Table 1: Time required for Cramer's rule
n 10 12 14 16 18 20
Time 0.4 sec. 1 min. 3.6 hrs. 41 days 38 years 16,000 years


Footnotes:

1Dahmen and Reusken, Numerik für Ingenieure und Naturwissenschaftler, Springer, 2006.


File translated from TEX by TTH, version 4.03.
On 15 Jul 2025, 22:55.