#07 (09/13/2023)

Newton-Raphson method (single variable)

The formula for the Newton-Raphson method can be derived by retaining the first two terms of the Taylor series of f(x) as
f(x) = f(xo) + f′(xo) (x − xo) + (higher order terms).
(1)
With f(x) = 0, the above equation can be solved for x as
x = xo − f(xo)
f′(xo)
 ,
(2)
which can be used as the iteration scheme:
For example, √2 can be computed by solving f(x) ≡ x2 − 2 = 0.
f(x) = x2 −2,     f′(x) = 2 x.
Start with x1 = 2.0 (initial guess), then,
x1
=
2,
x2
=
2− f(2)
f ′(2)
 = 2 − 2
4
 = 1.5,
x3
=
1.5 − f(1.5)
f ′(1.5)
 = 1.5 − 0.25
3
 = 1.41667,
x4
=
1.4167− f(1.4167)
f ′(1.4167)
 = … = 1.4142.
The convergence is attained after only 4 iterations.
More examples:
  1. Square root of a ( = √a)
    Let f(x) = x2 − a. It follows
    xn+1
    =
    xn − xn2 − a
    2 xn
    =
    1
    2
    		 ⎛
    ⎝     		xn + a
    xn
    		⎞
    ⎠     		.
    		(3)
    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/x − a. It follows
    xn+1
    =
    xn −
    1
    xn
    		 −a
    −1
    xn2
    =
    xn (2 − a xn).
    		(4)
    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/x2 − a, it follows

    xn+1
    =
    xn − 1/xn2 − a
    − 2/xn3
    =
    xn (3 − a xn2)/2.
    		(5)

Newton-Raphson method (multiple variables)

We want to solve a set of (non-linear) simultaneous equations in the form of
f1(x1, x2, x3, …xn)
=
0,
f2(x1, x2, x3, …xn)
=
0,
…
…
fn(x1, x2, x3, …xn)
=
0.
By expanding each of the above by the Taylor series, one obtains
f1(x)  ∼ f1 (xo) + ∂f1
∂x1
 ⎢
⎢
x0 
 (x1 − x10) + ∂f1
∂x2
 ⎢
⎢
x0 
 (x2 − x20) + …+ ∂f1
∂xn
 ⎢
⎢
x0 
 (xn − xn0),
f2(x)  ∼ f2 (xo) + ∂f2
∂x1
 ⎢
⎢
x0 
 (x1 − x10) + ∂f2
∂x2
 ⎢
⎢
x0 
 (x2 − x20) + …+ ∂f2
∂xn
 ⎢
⎢
x0 
 (xn − xn0),
…
fn(x)  ∼ fn (xo) + ∂fn
∂x1
 ⎢
⎢
x0 
 (x1 − x10) + ∂fn
∂x2
 ⎢
⎢
x0 
 (x2 − x20) + …+ ∂fn
∂xn
 ⎢
⎢
x0 
 (xn − xn0).
If x satisfies
f1(x) = 0,     f2(x) = 0,     …fn(x) = 0,
the above equations can be written as
⎛
⎜
⎜
⎜
⎜
⎜
⎝
0
0
:
0
⎞
⎟
⎟
⎟
⎟
⎟
⎠ 		= ⎛
⎜
⎜
⎜
⎜
⎜
⎝
f1
f2
:
fn
⎞
⎟
⎟
⎟
⎟
⎟
⎠ 		+ ⎛
⎜
⎜
⎜
⎜
⎜
⎜
⎜
⎜
⎜
⎝
∂f1
∂x1
 ,
∂f1
∂x2
 ,
…
∂f1
∂xn
∂f2
∂x1
 ,
∂f2
∂x2
 ,
…
∂f2
∂xn
…
…
…
…
∂fn
∂x1
 ,
∂fn
∂x2
 ,
…
∂fn
∂xn
⎞
⎟
⎟
⎟
⎟
⎟
⎟
⎟
⎟
⎟
⎠ 		⎛
⎜
⎜
⎜
⎜
⎜
⎝
x1 − x10
x2 − x20
:
xn − xn0
⎞
⎟
⎟
⎟
⎟
⎟
⎠ 		,
(6)
or

0 = f(xo) + J (x − xo),
(7)
where J is the Jacobian matrix defined as

J
≡
⎛
⎜
⎜
⎜
⎜
⎜
⎜
⎜
⎜
⎜
⎝
∂f1
∂x1
∂f1
∂x2
…
∂f1
∂xn
∂f2
∂x1
∂f2
∂x2
…
∂f2
∂xn
…
…
…
…
∂fn
∂x1
∂fn
∂x2
…
∂fn
∂xn
⎞
⎟
⎟
⎟
⎟
⎟
⎟
⎟
⎟
⎟
⎠
=
∂(f1, f2, …, fn)
∂(x1, x2, …, xn)
 .
(8)
Equation (7) can be solved for x as

x = xo − J−1 f(xo),
(9)
where J−1 is the inverse matrix of J.
Example:
Solve numerically for
x3 + y2 = 1,     x y = 1
2
 .
(10)
(Solution) Let
f1 ≡ x3 + y2 − 1,     f2 ≡ x y − 1
2
 .
It follows
J = ∂(f1, f2)
∂(x, y)
 = ⎛
⎜
⎜
⎝
3 x2,
2 y
y,
x
⎞
⎟
⎟
⎠ 		,
and
J−1 = ⎛
⎜
⎜
⎜
⎜
⎜
⎝
x
3 x3 − 2 y2
 ,
− 2 y
3 x3 − 2 y2
− y
3 x3 − 2 y2
 ,
3x2
3 x3 − 2 y2
⎞
⎟
⎟
⎟
⎟
⎟
⎠ 		,
so
J−1 f = ⎛
⎜
⎜
⎜
⎜
⎜
⎝
x
3 x3 − 2 y2
 ,
− 2 y
3 x3 − 2 y2
− y
3 x3 − 2 y2
 ,
3x2
3 x3 − 2 y2
⎞
⎟
⎟
⎟
⎟
⎟
⎠ 		⎛
⎜
⎜
⎝
x3 + y2 − 1
x y − 1/2
⎞
⎟
⎟
⎠ 		= ⎛
⎜
⎜
⎜
⎜
⎜
⎝
x4−x (y2+1)+y
3 x3−2 y2
4 x3 y−3 x2−2 y3+2 y
6 x3−4 y2
⎞
⎟
⎟
⎟
⎟
⎟
⎠ 		.
(11)
Hence the iteration scheme is
⎛
⎜
⎜
⎝
xn+1
yn+1
⎞
⎟
⎟
⎠ 		= ⎛
⎜
⎜
⎝
xn
yn
⎞
⎟
⎟
⎠ 		−⎛
⎜
⎜
⎜
⎜
⎜
⎝
xn4−xn (yn2+1)+yn
3 xn3−2 yn2
4 xn3 yn−3 xn2−2 yn3+2 yn
6 xn3−4 yn2
⎞
⎟
⎟
⎟
⎟
⎟
⎠ 		.
(12)
Sample C code: 
#include <stdio.h>
#include <math.h>
int main()
{
double x=1.0, y=1.0;
int i;

for (i=0;i<=10;i++)
  { x = x -(pow(x,4) + y - x*(1 + y*y))/(3*pow(x,3) - 2*y*y);
    y = y -(-3*x*x + 2*y + 4*pow(x,3)*y - 2*pow(y,3))/(6*pow(x,3) - 4*y*y);
  }

printf("%f %f\n", x, y);
return 0;
}
Sample Matlab/Octave code: 
x=1; y=1;
for i=0:10
 eqs=[ x*x*x+y*y-1; x*y - 1/2];
 jacob=[3*x*x 2*y; y x];
 right=inv(jacob)*eqs;
 x=x-right(1);
 y=y-right(2);
end;

fprintf('%f %f\n', x, y);
Sample Mathematica code: 
FindRoot[{x^3 + y^2 == 1, x y == 1/2}, {{x, 1}, {y, 1}}]
Online C compiler
Online Matlab/Octave
Wolfram Alpha
Starting (x, y) = (1.0, 1.0), the convergence was reached at (x, y) = (0.877275, 0.569947) after only 4 iterations. Note that this is just one of multiple roots. It is necessary to try different initial guess to obtain other roots.

Another example of Jacobian

Evaluate
I ≡ ⌠
⌡ 		1/2

0 
 ⌠
⌡ 		1 −2 y

0 
 exp⎛
⎝ 		x
x + 2 y
⎞
⎠ 		d x d y.
(13)
Hint:
  1. Draw the integral area.
  2. Change the variables from (x, y) to (u, v) by u = x, v = x + 2y.
    d x d y = ⎢
    ⎢     		∂(x, y)
    ∂(u, v)
    		⎢
    ⎢     		d u d v.
    		(14)

Implicit functions

If the relationship between x and y is given implicitly as

f(x, y) = 0,
(15)
it is still possible to express y as a function of x by Taylor series about x = a as

y(x) = y(a) + y′(a)(x − a) + y"(a)
2!
 (x − a)2 + y"′(a)
3!
 (x − a)3 + …
(16)
if y(a), y′(a), y"(a) … are known. One can differentiate f(x, y) = 0 with respect to x to get y′(x) without explicitly solving f(x, y) = 0 for y(x) as

fx (x, y) + fy(x, y) y′(x) = 0,
(17)
or

y′(x) = − fx(x,y)
fy(x,y)
 .
(18)
Once y′(x) is obtained, y"(x) can be derived by differentiating the both sides of y′(x) with respect to x again 1 as

y"(x)
=
− (fxx + fxy y′)fy − fx (fyx + fyyy′)
fy2
=
−
⎛
⎝ 		fxx + fxy ( − fx
fy
 )⎞
⎠ 		fy − fx ⎛
⎝ 		fyx + fyy( − fx
fy
 )⎞
⎠
fy2
=
2fx fy fxy −fx2 fyy − fy2 fxx
fy3
 ,
(19)
and other higher order derivatives can be obtained in a similar manner.
It is seen that for y = y(x) to exist at x = x0, the denominator of fy must not vanish, i.e.
∂f
∂y
 (x0, y0) ≠ 0.
(20)

Example

Expand y(x) about x = 1 up to and including second order for
f(x, y) ≡ x y −y3 − 1 = 0.
(21)
Solution: The Taylor series of y(x) about x=1 is
y(x) = y(1) + y′(1)(x − 1) + y"(1)/2! (x − 1)2 + …
(22)
By substituting x = 1 into x y − y3 − 1 = 0, one gets y(1) = −1.32 … by solving y − y3 − 1 = 0 (there are two other complex roots but are not considered here). Next, by differentiation x y − y3 − 1 = 0 with respect to x (note that y is a function of x), one gets y + x y′− 3 y2 y′ = 0 which can be solved for y′ as y′ = y/(3 y2 − x). This can be further differentiated with respect to x as y" = (y′(3y2 − x) − y(6 y y′−1))/(3 y2 − x)2 so it is possible to obtain y(1), y′(1), y"(1), ... sequentially. Finally, one obtains

y(x) = −1.32472 −0.310629 (x−1) +0.0170796 (x−1)2+0.00114503 (x−1)3 −0.00128188 (x−1)4+ …
(23)
implicit_diff.jpg

Footnotes:

1
⎛
⎝ 		u
v
⎞
⎠ 		′

 
 = u′v − u v′
v2
 .


File translated from TEX by TTH, version 4.03.
On 12 Sep 2023, 21:27.