Modify newton from the previous exercise so that it works on

Modify newton from the previous exercise so that it works on a system of equations F(x). The function fprime now returns the Jacobian matrix J, and the Newton update is written mathematically as although in numerical practice one does not compute the inverse of the Jacobian but solves a linear system of equations in which F(xn) is the right-hand side.

Solution

error1 = 1.e8;
xx(1) = 0; % initial guesses
xx(2) = 0;
iter=0;
maxiter=30.
% begin iteration
while error1>1.e-12
iter=iter+1;
x = xx(1);
y = xx(2);
% calculate the functions
f(1) = 10*x+3*y*y-3;
f(2) = x*x-exp(y) -2.;
% calculate the Jacobian
J(1,1) = 10;
J(1,2) = 6*y;
J(2,1) = 2*x;
J(2,2) = -exp(y);
% solve the linear equations
y = -J\\f\';
% move the solution, xx(k+1) - xx(k), to xx(k+1)
xx = xx + y\';
% calculate norms
error1=sqrt(y(1)*y(1)+y(2)*y(2))
error(iter)=sqrt(f(1)*f(1)+f(2)*f(2));
ii(iter)=iter;
if (iter > itermax)
error1 = 0.;
s=sprintf(\'****Did not converge within %3.0f iterations.****\',itermax);
disp(s)
end
% check if error1 < 1.e-12
end
x = xx(1);
y = xx(2);
f(1) = 10*x+3*y*y-3;
f(2) = x*x-exp(y) -2.;
% print results
f
xx
iter
% plot results
semilogy(ii,error)
xlabel(\'iteration number\')
ylabel(\'norm of functions\')