kill(all);

load(draw);

r1(x,y) := x^2 - x + y-1/4;

r2(x,y) := exp(x+y)-y*x;

wxplot3d(r1(x,y),[x,-1,1],[y,-1,1]);

wxplot3d(r2(x,y),[x,-1,1],[y,-1,1]);

wxcontour_plot([r1(x,y),r2(x,y),[x,-1,1],[y,-1,1]],[gnuplot_preamble, "set cntrparam levels incremental -4,1,4"]);

definizione di funzione vettore di variabile vettore

r(u) := [r1(u[1],u[2]),r2(u[1],u[2])];

definizione dello jacobiano

define(J(u), 
    (
        tmp : jacobian([r1(x,y),r2(x,y)],[x,y]),
        ev( tmp , x=u[1],y=u[2] )
    )
);

definizione dello jacobiano inverso

define(invJ(u), 
    (
        tmp : jacobian([r1(x,y),r2(x,y)],[x,y]),
        tmp : invert(tmp),
        ev( tmp , x=u[1],y=u[2] )
    )
);

definizione del passo N-R

NRnext(u) := ev(u - list_matrix_entries(invJ(u).r(u)),numer);

definisco tre iterati N-R a partire da un u0 di partenza

u0 : [0.5,-0.3];

u1 : NRnext(u0);
u2 : NRnext(u1);
u3 : NRnext(u2);

cp1 : [
color=red,
line_width=2,
implicit(r1(x,y)-0,x,-1,1,y,-1,1),
line_width=1,
implicit(r1(x,y)-1,x,-1,1,y,-1,1),
implicit(r1(x,y)-2,x,-1,1,y,-1,1),
color=dark_red,
implicit(r1(x,y)+1,x,-1,1,y,-1,1)
]$
cp2 : [
color = green,
line_width=2,
implicit(r2(x,y)-0,x,-1,1,y,-1,1),
line_width=1,
implicit(r2(x,y)-4,x,-1,1,y,-1,1),
implicit(r2(x,y)-3,x,-1,1,y,-1,1),
implicit(r2(x,y)-2,x,-1,1,y,-1,1),
implicit(r2(x,y)-1,x,-1,1,y,-1,1)
]$
vecs :  [color=blue,vector(u0,u1-u0),vector(u1,u2-u1),vector(u2,u3-u2)]$

wxdraw2d(cp1,cp2,head_length=0.05,vecs);

wxcontour_plot([r1(x,y),r2(x,y),[x,-1,1],[y,-1,1]],[gnuplot_preamble, "set cntrparam levels incremental -4,1,4"]);