Implementation of Newton's method
for a general system of equations
In this exercise you will implement Newton's method. This program
will be used in Project 12.
1.
Reading:
90. Linearization and Newton's method.
(ps)
The teacher presents:
the derivative (Jacobi matrix) Df(x) of a function f:
Rn ->Rm, numerical computation of the derivative,
Newton's method.
2. Write a program that computes the Jacobi matrix of a function according to the template jacobi.m .
Consider first a function of one variable f: R ->R. Use the central difference quotient: ( f(x+h ) - f(x-h ) ) / (2h) to approximate the derivative with respect to x .
The matlab command feval(f,x) evaluates the function
whose name is stored in the variable f at the point x.
Examples:
>> feval('sin',pi)
>> funk='sin', x=pi/2, h=.1
>> feval(funk,x)
>> feval(funk,x+h)-feval(funk,x-h)
Try this on several examples until you
understand it!! Don't continue until you really understand
how this command works.
Debugging (finding mistakes): remove all semicolons (;) in the program. Choose a simple example and compute everything by hand and compare with what the program computes. Use this technique all the time! Compute with pen and paper all the time to check what the computer does.
Test your program!
>> z = jacobi('sin', 0)
>> z = jacobi('sin', pi/2)
>> z = jacobi('sin', pi)
Now consider a function of several variables f: Rn ->Rm. Use the central difference quotient: ( f(x+h ei ) - f(x-h ei ) ) / (2h) to approximate the partial derivative with respect to xi . Here ei is the unit vector in the ith direction.
This formula computes the ith column of the Jacobi matrix.
First create a zero vector e of the same size as x. Then change the the ith element of e to 1. Hint: zeros, size.
In matlab A(:,i) is column number i of the matrix A. Try this on some matrices until you understand it!!
Test your Jacobi program on several functions f and several points
x:
compute the derivative by hand and then run the program. Don't stop
until you are sure that the program works correctly . What
happens if you type
>> help jacobi ?
How large is
the error in each case? How large should the step h
be chosen? See AMBS 23.13.
Examples: (in matlab notation)
sin(x), x=0, x=pi/2, x=pi
f(x)=x'*x, x=[1;2;3]
f(t)=[t;t^2;t^3], t=0, t=1
f(x)=x, x=[7;-9;3;5]
f(x)=x'*x, x=[7;-9;3;5]
f(x)=[ x(2)*(1-x(1)^2); 2-x(1)*x(2) ], x=[1;2]
f(x)=[ x(2)*(1-x(1)^2)+x(3)^2; 2-x(1)*x(2) ], x=[1;2;3]
You find more examples among the problems in 90. Linearization and Newton's method.
3. Write a program that implements newton's method according to the template newton.m .
while |h|< tol
evaluate the Jacobian:
A=Df(x)
evaluate the residual:
b=-f(x)
solve the linearized equation: Ah=b
(use matlab's slash command: h=A\b )
update:
x=x+h
Try each command separately until you understand it!! Don't put it into the newton program before you are sure how it works.
Debugging: remove all semicolons (;) in the program, set the tolerance very large (tol=10) so that the program will do only one iteration. Choose a simple example and compute everything by hand and compare with what the program computes.
Test your program on several examples: solve the equations by hand, then run the program. Make sure that you find all solutions.
Begin with the following examples. Then find some test equations
of your own.
f(x)=sin(x)
f(x)=x^2-2
f(x)=[ x(2)*(1-x(1)^2); 2-x(1)*x(2) ]
f(x)=[ x(1)*(1-x(2)); x(2)*(1-x(1)) ]
f(x)=[ x(1)+x(1)*x(2)^2+x(1)*x(3)^2; -x(1)+x(2)-x(2)*x(3)+x(1)*x(2)*x(3);
x(2)+x(3)-x(1)^2 ]
See also the problems in 90. Linearization and Newton's method.
/stig