On 02/14/2013 01:13 PM, Hugo Coolens wrote:
> I'm looking for the extrema of a function w=f(x,y,z). Can someone here
> point me to how to accomplish this with Maxima?
>
> thanks in advance
> hugo
>
Hi Hugo,
since I like dynamical systems, I would interpret your problem as the
problem of finding the equilibrium points of a dynamical system and
characterizing them.
Here is a specific example:
(%i1) w: 2*exp(-(x-3)^2-(y-1)^2-(z-1)^2) + exp(-(x-1)^2-(y-5)^2-(z-2)^2)$
define a "gradient system"
(%i2) g: [diff(w,x),diff(w,y),diff(w,z)]$
g is the phase velocity of that system in the (x,y,z) phase space. The
critical points of w will now be the equilibrium points of that
dynamical system; thus, all the points that solve g=0 must be found.
That is a difficult task which, in general, can only be done numerically.
(%i3) load("mnewton");
in this case it is easy to see that there must be 3 equilibrium points
near (3,1,1), (1,5,2), and somewhere near the middle point of the
segment p1p2, (2,3,1.5)
(%i4) p1: first(mnewton(g,[x,y,z],[3,1,1]));
(%o4) [x = 2.999999999241744, y = 1.000000001516512,
z = 1.000000000379128]
(%i5) p2: first(mnewton(g,[x,y,z],[1,5,2]));
(%o5) [x = 1.000000003033024, y = 4.999999993933951,
z = 1.999999998483488]
(%i6) p3: first(mnewton(g,[x,y,z],[2,3,1.5]));
(%o6) [x = 1.963516864231185, y = 3.07296627153763,
z = 1.518241567884407]
To characterize each of those 3 equilibrium points, we now calculate the
jacobian matrix for the system (which is also the hessian matrix of w)
(%i7) H: jacobian(g,[x,y,z])$
since we are dealing with a gradient system, its jacobian matrix is
symmetric, so all of its eigenvalues will be real (either nodes or
saddle points for the dynamical system, but no centers or foci). The
eigenvalues of a symmetric matrix can be calculated using eigens_by_jacobi
(%i8) H1: subst(p1,H)$
(%i9) eigens_by_jacobi(H1);
(%o9) [[- 4.000000001516512, - 3.999999937823004,
- 4.000000001516512], ...
since the three eigenvalues are all negative, then p1 is an attractive
node (a local maximum of w).
(%i11) eigens_by_jacobi(H2);
(%o11) [[- 2.000000003033024, - 1.999999875646001,
- 2.000000003033024], ...
p2 is also a local maximum
(%i14) eigens_by_jacobi(H3);
(%o14) [[- .02949737665668998, .2798128317058039,
- .02949737665668998], ....
p3 is a saddle point.
I hope this helps you.
Best regards,
Jaime