CODEWITHSUNDEEP
mathevaluationmaximanewtons-method
useslessone
useslessone

Reputation: 69

maxima: evaluating Jacobian at a point

I'm not very good with maxima, however I'm trying to learn.
I am attempting to write a Newton iterative solver for a fairly small nonlinear system. In order to do this I must evaluate the Jacobian at the current iteration. But, I can't seem to figure out a way to evaluate the Jacobian at all.

I currently have:

F1000(D5,D6,D7,D8,D9,D10,b1,b2,b3,b4,b5,b6,b7,b8,b9,d1,d2,d3,e1,e2,e3,f1,f2,f3) 
     := (-f1-3*d1-9*b1+60)*l1$
F1001(D5,D6,D7,D8,D9,D10,b1,b2,b3,b4,b5,b6,b7,b8,b9,d1,d2,d3,e1,e2,e3,f1,f2,f3) 
     := (-f2+e2-6*d2-3*b2-15)*l1$
/*further functions omitted*/

J : jacobian([ 
   F1000(D5,D6,D7,D8,D9,D10,b1,b2,b3,b4,b5,b6,b7,b8,b9,d1,d2,d3,e1,e2,e3,f1,f2,f3),
   F1001(D5,D6,D7,D8,D9,D10,b1,b2,b3,b4,b5,b6,b7,b8,b9,d1,d2,d3,e1,e2,e3,f1,f2,f3),
   /*functions omitted*/
   F0018(D5,D6,D7,D8,D9,D10,b1,b2,b3,b4,b5,b6,b7,b8,b9,d1,d2,d3,e1,e2,e3,f1,f2,f3),
   F0019(D5,D6,D7,D8,D9,D10,b1,b2,b3,b4,b5,b6,b7,b8,b9,d1,d2,d3,e1,e2,e3,f1,f2,f3)
   ],
   [D5,D6,D7,D8,D9,D10,b1,b2,b3,b4,b5,b6,b7,b8,b9,d1,d2,d3,e1,e2,e3,f1,f2,f3]
   )$

But, now i cannot find a good way to evaluate J at, for example, 

u_init : [D5_0,D6_0,D7_0,D9_0,...,e3_0,f1_0,f2_0,f3_0]$

I've tried things like:

subs(J, [D5,D6,D7,D8,D9,D10,b1,b2,b3,b4,b5,b6,b7,b8,b9,d1,d2,d3,e1,e2,e3,f1,f2,f3],
     [1,1,1,1...,1,1,1,1]),

and 

ev(J, [D5,D6,D7,D8,D9,D10,b1,b2,b3,b4,b5,b6,b7,b8,b9,d1,d2,d3,e1,e2,e3,f1,f2,f3],
     [1,1,1,1...,1,1,1,1]),

and I just cannot find a good way to evaluate the Jacobian, $J$, at the point $u_i$. I really appreciate any clarifications/insights.

Thanks

Upvotes: 2

Views: 3443

Answers (1)

Roberto
Roberto

Reputation: 3083

Try

subst([D5=D5_0,D6=D6_0,...,f3=f3_0], J);

An example for a smaller Jacobian:

(%i1) F1(x,y):=x^2+y;
                                           2
(%o1)                         F1(x, y) := x  + y
(%i2) F2(x,y):=-x-y^2;
                                                2
(%o2)                        F2(x, y) := - x - y
(%i3) J:jacobian([F1(x,y),F2(x,y)],[x,y]);
                                [ 2 x    1   ]
(%o3)                           [            ]
                                [ - 1  - 2 y ]
(%i4) subst([x=2,y=3],J);
                                 [  4    1  ]
(%o4)                            [          ]
                                 [ - 1  - 6 ]

Upvotes: 2

Related Questions