Problem concerning Numerical Solutions of Nonlinear Systems of Equations (Burden and Faires)

Question

I need help with this particular question:

The nonlinear system

1. $3x_1 - \cos (x_2 x_3) - \frac{1}{2} = 0$
2. $x{_1}^2 - 625x{_2}^2 - \frac{1}{4}=0$
3. $\exp ^{-x_1x_2} + 20x_3 + \frac{10\Pi -3}{3}=0$

has a singular Jacobian matrix at the solution. Apply Newton's method with x$^{(0)} = (1,1,-1)^t$. Note that convergence may be slow or may not occur within a reasonable number of iterations.

My attempt on Maple:

What am I doing wrong here?

Solving it by-hand is a pain.

Answer 1

In your NewtonMV procedure, try changing,

 A := subs(s,f);

to,

 A := eval(J, [s]);

or to,

 A := subs(s, J);

Also, change the call LinearSolve(...) to :-LinearAlgebra:-LinearSolve(...) so that it gets the right export from the right top package, regardless of whether you've loaded any subpackage of Student.

Personally, I'd rather see evalf wrapped around each of the two arguments to LinearSolve than around its call, if only so that the very first iteration doesn't attempt an exact subcomputation.

So, something like,

 dx := :-LinearAlgebra:-LinearSolve( evalf(A), evalf(-fx) );

Please note that you can inline code into your Questions and Answers here see the menubar above the editing region. It's not so helpful to paste images of code.

[edited, after all the typos have been resolved] You might want to consider additional stopping criteria. Apart from a maximal number of iterations to allow, you could stop when either the increment dx got small enough, or when the evaluation of f got close enough to zero (ie. x is considered a root).

restart:

NewtonMV:=proc(f,v,x0,n,xtol,ftol)
local i,result,J,x,fx,A,dx,s;
J:=VectorCalculus:-Jacobian(f,v);
x:=x0;
dx:=Vector([infinity,infinity,infinity]);
fx:=Vector([infinity,infinity,infinity]);
result:=convert(x0,list);
for i from 1 to n
  while :-LinearAlgebra:-Norm(dx)>xtol
    and :-LinearAlgebra:-Norm(fx)>ftol do
   s:=[seq(v[i]=x[i],i=1..nops(v))];
   fx:=evalf(eval(f,s));
   A:=evalf(eval(J,s));
   dx:=:-LinearAlgebra:-LinearSolve(evalf(A),evalf(-fx));
   if indets(dx,name)<>{} then
      dx:=subs([seq(p=0,p=indets(dx,name))],dx);
   end if;
   x:=x+dx;
   result:=result,convert(x,list);
end do;
printf("\nnumber of iterations: %a\n\n",i-1);
printf("final increment dx, and its norm:\n      %a ,  %a\n\n",
       convert(dx,list), :-LinearAlgebra:-Norm(dx));
printf("substitution of solution x into f, and its norm:\n      %a ,  %a\n\n",
       convert(eval(f,s),list), :-LinearAlgebra:-Norm(fx));
return [result];
end proc:

lhs1:=3*x1-cos(x2*x3)-1/2;
lhs2:=x1^2-625*x2^2-1/4;
lhs3:=exp(-x1*x2)+20*x3+(10*Pi-3)/3;

NewtonMV(Vector([lhs1,lhs2,lhs3]),[x1,x2,x3],
         Vector([1,1,-1]), 400, 1e-5, 1e-7);

That is not intended to be the epitome of great coding, but it might give you some ideas.