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.
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.
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