) > 0
UCES
Course: Methods and Analysis in UCES
Applied Area: falling mass, heat transfer with heat source
Method: Newton's algorithm, Euler-Trapezoid method with Newton solver
Assessment: quadratic and local convergence theorem
Prerequisites:
Objectives:
General Information: We establish quadratic and local convergence of Newton's algorithm for a single function of one variable. This will later be generalized to N functions and N unknowns.
Contact Person: R. E. White, NCSU, white@math.ncsu.edu
Revision Date: 8-3-94
Copyright ©1994, R. E. White. All rights reserved.
Newton's algorithm is one of the most important numerical schemes because, under appropriate conditions, it has local convergence and quadratic convergence properties. Local convergence means that if the initial guess is sufficiently close to the root, then the algorithm will converge to the solution. This is in contrast to the Picard algorithm (recall our attempts to find the square root of 2). Quadratic convergence means that the error at the next step will be proportional to the square of the error at the current step. In general, the Picard algorithm only has first order convergence where the error at the next step is proportional to just the error at the present step (see the proof of the Picard convergence theorem).
We give two applications. The Euler-Trapezoid algorithm is used to solve the stiff differential equation which models a mass falling through a viscous medium. A second application evolves from continuum solution to a Newton's law of cooling model with heat source.
Application to a Falling Mass Through a Viscous Medium. Consider a mass which is falling through a viscous medium. There will be two forces: one from gravity, mg, and a second from the viscous medium, -ku2. Here the speed of the mass is denoted by u, and the k is a constant which reflects the mass and the viscosity of the medium. Newton's law of motion then gives the differential equation
If k or u is initially very large, then this will be a stiff differential
equation. As time evolves, the speed will approach
In the Euler-Trapezoid algorithm we can use Newton's algorithm in place of the Picard algorithm which may require some constraint on the step size. We must replace the fixed point problem u = g(u) by a root problem F(u) = g(u) - u = 0. The inner loop of the Euler-Trapezoid algorithm should be replaced by the Newton scheme to solve
For the initial guess, we can choose the previous speed u(k) which should be relatively close to this root.
Application to Newton's Law of Cooling with Heat
Source. Consider the differential equation
This has at least one solution because F(0) = 1 - e < 0 and
) > 0
The model has either the form of a fixed point x = g(x) or a root of f(x) = 0 where f or g are given functions. Newton's method is formulated for the root problem, but clearly the fixed point problem can be written as a root problem when one defines f(x) = x - g(x).
The idea behind Newton's algorithm is to approximate the function f(x) at a given point by a straight line. Then find the root of the equation associated with this straight line. One continues to repeat this until little change in approximated roots is observed.
 |
| Figure: Newton's Algorithm |
The equation for the straight line at the mth iteration is
Define xm+1 so that y = 0 where the straight line intersects the x axis
| xm+1 = xm - f(xm )/f'(xm ). | (1) |
There are two common stopping criteria for Newton's algorithm. First, is the common absolute error which requires two successive iterates to be close. The second measures how close the function is to zero.
choose xold, eps1 > 0 and eps2 > 0
for m = 1, maxit
xnew = xold - f(xold)/f'(xold)
if |xnew - xold| < eps1 and |f(xnew)| < eps2 go to end
xold = xnew
endloop
end.
Newton's Algorithm for Euler-Trapezoid Method.
for m = 1, maxit
F = u(k) + (h/2)*(oldf + f((k+1)*h,oldu)) - oldu
FP = (h/2) fu((k+1)*h,oldu) - 1
newu = oldu - F/FP
if |oldu - newu| < eps go to end
oldu = newu
endloop.
Example. Consider f(x) = 2 - x2 = 0. The derivative of f is
-2x, and the iteration in (1) can be viewed as a special Picard algorithm where
| m | xm | error = Em | Em/(E(m-1)2 ) |
| 0 | 1.0 | .414213 | |
| 1 | 1.5 | .085786 | 2.000005 |
| 2 | 1.4166666 | .002453 | 3.000097 |
| 3 | 1.4142156 | .000002 | 3.008604 |
Example. Consider f(x) = x1/3 - 1. The derivative of f is
(1/3)x-2/3. The corresponding Picard iterative function is
| x0 = initial guess | m = iter. for conv. | x0 = initial guess | m = iter. for conv. | |
| 10.0 | no conv. | 00.1 | 4 | |
| 05.0 | no conv. | -0.5 | 6 | |
| 04.0 | 6 | -0.8 | 8 | |
| 03.0 | 5 | -0.9 | 20 | |
| 01.8 | 3 | -1.0 | no conv. |
We consider the differential equation for the mass falling through a viscous
medium which exerts a force -.5u2 where u is the speed of the mass.
If the mass is one unit, then the equation of motion is
Maple Code for Euler-Trapezoid Algorithm with Newton Solver.
Input Data:
> t:='t':
> u:='u':
> f:=(t,u)->32 -.5*u*u;
> fp:=(t,u)->-u;
> t0:=0:
> u0:=100:
> T:=.20:
> K:=10:
Define Procedure EULERTN:
> EULERTN:=proc(f,fp,t0,u0,T,K)
> t(0):= t0:
> h:=T/K:
> u(0):= u0:
> eps:=.0001:
> maxm:= 50:
> for k from 0 to K do
> oldf:=f(t(k),u(k)):
> oldu:=u(k):
> t(k+1):= t(k) + h:
> test:=1:
> for m from 1 to maxm while(test > eps) do
> F:=u(k) + (h/2)*(oldf + f(t(k+1),oldu)) - oldu:
> FP:= (h/2)*fp(t(k+1),oldu)-1:
> newu:= oldu - F/FP:
> test:= abs(oldu - newu):
> oldu:=newu:
> od:
> u(k+1) :=newu:
> od:
> end:
Output Data:
> with(plots):
> EULERTN(f,fp,0,100,.2,10):
> L:=[seq([t(n),u(n)],n=0..10)]:
> EULERTN(f,fp,0,100,.20,20):
> LL:=[seq([t(n),u(n)],n=0..20)]:
> EULERTN(f,fp,0,100,.20,40):
> LLL:=[seq([t(n),u(n)],n=0..40)]:
> plot({L,LL,LLL},x=0..(.2));
 |
| Figure: Speed versus Time: K = 10, 20, 40 |
In the above examples we noted that Newton's algorithm was a special case
of the Picard algorithm with
| g'(x) = 1 - ((f'(x))2 - f(x)f''(x))/(f'(x)2 ) = f(x)(f''(x)/(f'(x)2 ). | (2) |
If x is a solution of f(x) = 0 and f(x) is continuous, then we can make f(x) as small as we wish by choosing x close to x. So, if in (2) f''(x)/(f'(x)2 ) is bounded, then g(x) will be contractive for x near x. Under the conditions listed in the following theorem this establishes the local convergence.
Newton's Convergence Theorem. Consider f(x) = 0 and assume x is such a root. If f'(x) is not zero, f''(x) is continuous on an interval containing x, then
| xm+1 - x | 
[max |f''(x)|/(2 min |f'(x)|)] |xm -
x|2. |
(3) |
Proof. In order to prove the quadratic convergence in (3), we use the
extended mean value theorem where
Divide the above by f'(xm ) and use the definition of xm+1 in (1) to get
| 0 | = -(xm+1 - xm) + (x - xm) + (f''(c)/(2f'((xm)) (x - xm)2 |
| = -(xm+1 - x) + (f''(c)/(2f'((xm)) (x - xm)2. |
Since f'(x) for some interval about x must be bounded away from zero, and f''(x) and f'(x) are continuous, the inequality in (3) must hold.