4

I have 3D Gaussians and their derivatives (Laplace operator)

f1(x,y,z):=exp(-w1*((x-x1)^2+y^2+z^2));
f2(x,y,z):=exp(-w2*(x^2+y^2+z^2));
dx2_f2(x,y,z):=diff(f2(x,y,z),x,2);
Lf2(x,y,z):=diff(f2(x,y,z),x,2) + diff(f2(x,y,z),y,2) + diff(f2(x,y,z),z,2);

Now I want to plot a profile along x cutting the function for fixed y,z for some values of parameters

w1:1.2;w2:0.5;x1:1.5;
plot2d( Lf2(x,0,0) ,[x,-5,5]);

I get following error

diff: second argument must be a variable; found 0
#0: Lf2(x=x,y=0,z=0)
 -- an error. To debug this try: debugmode(true);

While plot2d( [f1(x,0,0),f2(x,0,0),dx2_f2(x,0,0)] ,[x,-5,5]); works just fine.

The same error is if I try to partially substite the variables manually Lf2x(x):=Lf2(x,0,0);

My guess, the problem is Maxima/lisp does some lazy evaluation, therefore the derivatives along y,z are not yet calculated when I substitute them by y=0,z=0. Therefore it refuse do derivative along constant (?).

But I don't know how to solve it (i.e. substitute the constant only after the derivatives are calculated)

Prokop Hapala
  • 2,424
  • 2
  • 30
  • 59
  • 1
    The body of a function is not evaluated when the function is defined by ":=", so the derivatives aren't calculated until the function is called, when x, y, and z have numerical values. One way to solve this is to define the function like this: `define(dx2_f2(x,y,z), diff(f2(x,y,z),x,2))` because `define` evaluates the function body at the time the function is defined. Try `fundef(dx2_f2)` to see what Maxima thinks the function body is. – Robert Dodier May 29 '20 at 20:54
  • Thanks. And what about the `ev` and `''` operators? Can it help when function is already defined as `:=` ? I would prefer use `:=` since it looks more like math and it is more human readable. http://maxima.sourceforge.net/docs/manual/maxima_8.html – Prokop Hapala May 30 '20 at 07:11

1 Answers1

1

Here is a solution which makes use of quote-quote as you suggested. The main idea is to say foo(x) := ''(diff(something, x)) instead of foo(x) := diff(something, x).

(%i2) f1(x, y, z) := exp(-w1*((x - x1)^2 + y^2 + z^2));
                                           2    2    2
(%o2)   f1(x, y, z) := exp((- w1) ((x - x1)  + y  + z ))
(%i3) f2(x, y, z) := exp(-w2*(x^2 + y^2 + z^2));
                                        2    2    2
(%o3)       f2(x, y, z) := exp((- w2) (x  + y  + z ))
(%i4) dx2_f2(x, y, z) := ''(diff(f2(x, y, z), x, 2));
                                           2    2    2
                             2  2   - w2 (z  + y  + x )
(%o4) dx2_f2(x, y, z) := 4 w2  x  %e
                                                     2    2    2
                                              - w2 (z  + y  + x )
                                     - 2 w2 %e
(%i5) Lf2(x, y, z) := ''(diff(f2(x, y, z), x, 2) + diff(f2(x, y, z), y, 2) + diff(f2(x, y, z), z, 2));
                                        2    2    2
                          2  2   - w2 (z  + y  + x )
(%o5) Lf2(x, y, z) := 4 w2  z  %e
                     2    2    2
       2  2   - w2 (z  + y  + x )
 + 4 w2  y  %e
                     2    2    2                  2    2    2
       2  2   - w2 (z  + y  + x )          - w2 (z  + y  + x )
 + 4 w2  x  %e                    - 6 w2 %e
(%i6) w1: 1.2;
(%o6)                          1.2
(%i7) w2: 0.5;
(%o7)                          0.5
(%i8) x1: 1.5;
(%o8)                          1.5
(%i9) plot2d( Lf2(x, 0, 0) ,[x, -5, 5]);
(%i10) plot2d( [f1(x, 0, 0), f2(x, 0, 0), dx2_f2(x, 0, 0)], [x, -5, 5]);

With these definitions, I find the plots at the end make nice looking plots without errors. Note that the function definitions %o4 and %o5 have the actual derivatives on the right-hand side, not a diff expression. You could get the same effect by using define as I suggested at first.

Robert Dodier
  • 16,905
  • 2
  • 31
  • 48