Implementation of Newton's method for a general system of equations
1. 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 Lecture 1.
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 Lecture 1.
4. Write the documentation of your work for your individual folder. The documentation must contain: a printout of the program (well written and with comments), a page describing what the program does, what the algorithm is, a couple of good test examples. All this should be well written and informative, not too long, not too short. Take it seriously; I will read it when I grade the final exam.
Note: each student must keep an individual folder containing documentation of the obligatory work such as programs and projects. The folder must be handed in at the exam. Each student must also keep a copy of the each obligatory program that we write. These programs may be used in future courses.
Have fun!
/stig