The Lagrange Method
>f &= x^(4/5)*y^(1/5)
4/5 1/5
x y
>g &= x+y-1
y + x - 1
>L &= f+c*g
4/5 1/5
c (y + x - 1) + x y
>&solve(gradient(L,[x,y,c]),[x,y,c]); sol &= %[1]
8/5
4 1 2
[x = -, y = -, c = - ----]
5 5 5
>plot2d(f,niveau=f(4/5,1/5),a=0,b=2,c=0,d=2); ...
>plot2d(g,niveau=0,add=1,color=5); insimg;

>&f with sol
4/5
4
----
5
>H &= matrix([0,diff(g,x),diff(g,y)], ...
> [diff(g,x),diff(L,x,2),diff(diff(L,x),y)], ...
> [diff(g,y),diff(diff(L,x),y),diff(L,y,2)])
[ 0 1 1 ]
[ ]
[ 1/5 ]
[ 4 y 4 ]
[ 1 - ------- ------------ ]
[ 6/5 1/5 4/5 ]
[ 25 x 25 x y ]
[ ]
[ 4/5 ]
[ 4 4 x ]
[ 1 ------------ - ------- ]
[ 1/5 4/5 9/5 ]
[ 25 x y 25 y ]
>&H with sol, &:determinant(%)
[ 0 1 1 ]
[ ]
[ 4/5 ]
[ 1 4 ]
[ 1 - ------ ---- ]
[ 1/5 5 ]
[ 5 4 ]
[ ]
[ 4/5 9/5 ]
[ 4 4 ]
[ 1 ---- - ---- ]
[ 5 5 ]
3.78929141628
>&solve(g,y), &f with %, solx &= solve(diff(%,x)=0,x)
[y = 1 - x]
1/5 4/5
(1 - x) x
4
[x = -]
5
>&solve(g,y); &f with %; &diff(%,x,2) with solx | float
- 3.789291416275997
>&solve(g,y); &f with %; &diff(%,x,2)|factor
4
--------------------------
4/5 6/5
25 (1 - x) (x - 1) x
>Lag &= f+c*g with c=(c with sol), &eigenvalues(hessian(Lag,[x,y]) with sol)
8/5
4/5 1/5 2 (y + x - 1)
x y - ----------------
5
17
[[- ------, 0], [1, 1]]
1/5
5 4
>plot3d("Lag(x,y)",r=0.01,cx=0.8,cy=0.2, ...
> angle=160°,height=50°,>contour); insimg;

>function lg([x,y,c]) ...
x=abs(x); y=abs(y);
return [&:"diff(@f+c*@g,x)", ..
&:"diff(@f+c*@g,y)", ..
&:"diff(@f+c*@g,c)"]
endfunction
>broyden("lg",[1,1,-1])
[ 0.8 0.2 -0.606286626604 ]
>function h(v) := norm(lg(v))
>neldermin("h",[1,1,0])
[ 0.800000000001 0.2 -0.606286626604 ]
>&solvelagrange(f,g,[x,y])
8/5
2 1 4
[[lambda = ----, y = -, x = -]]
5 5 5
>sol &= solvelagrange(x*y^2*z^3,x+y+z-1,[x,y,z])
[[x = %r1, lambda = 0, z = 0, y = 1 - %r1],
[x = %r2, lambda = 0, z = 1 - %r2, y = 0],
1 1 1 1
[x = -, lambda = --, z = -, y = -],
6 72 2 3
1 2
[x = -, lambda = 0, z = 0, y = -], [x = 0, lambda = 0, z = 1, y = 0]]
3 3
>&[x,y,z] with sol[3]
1 1 1
[-, -, -]
6 3 2
>&getlagrange(x*y^2*z^3,x+y+z-1,[x,y,z])
2 3 3 2 2
[y z = lambda, 2 x y z = lambda, 3 x y z = lambda,
z + y + x - 1 = 0]
>function g(x,y) &= x^2*y+y^2*x+x^4-1
2 2 4
x y + x y + x - 1
>plot2d("g",niveau=0,r=2); insimg;

>function f(x,y) &= (1-x)^2+(1-y)^2
2 2
(1 - y) + (1 - x)
>L &= f(x,y)+lambda*g(x,y)
2 2 4 2 2
(x y + x y + x - 1) lambda + (1 - y) + (1 - x)
>function dL([x,y,lambda]) &= grad(L)
2 3
[(y + 2 x y + 4 x ) lambda - 2 (1 - x),
2 2 2 4
(2 x y + x ) lambda - 2 (1 - y), x y + x y + x - 1]
>function DdL([x,y,lambda]) &= jacobian(dL(x,y,lambda),[x,y,lambda])
[ 2 2 3 ]
[ (2 y + 12 x ) lambda + 2 (2 y + 2 x) lambda y + 2 x y + 4 x ]
[ ]
[ 2 ]
[ (2 y + 2 x) lambda 2 x lambda + 2 2 x y + x ]
[ ]
[ 2 3 2 ]
[ y + 2 x y + 4 x 2 x y + x 0 ]
>newton2("dL","DdL",[1,1,0]), plot2d(%[1],%[2],points=true,add=true); insimg;
[ 0.661468492772 0.823282064697 0.23150454352 ]

>&solve(dL(x,y,lambda))[6]
[lambda = 0.23150453851715, y = 0.82328207592504,
x = 0.66146848602989]
>{xsol,valid}=inewton2("dL","DdL",[1,1,0]); xsol,
[ ~0.6614684927715179,0.6614684927715191~
~0.8232820646967926,0.8232820646967943~
~0.2315045435202875,0.2315045435202889~ ]
>valid,
1
Examples Homepage