偏微分の変数変換

これも荒いが、下書きが増えすぎているので公開してしまおう。

EMAN 氏の http://homepage2.nifty.com/eman/analytic/bibun.html を Maxima でなぞってみる。なお、実用的には以前の記事の方がはるかに簡単。

define(x(r,t,p), r*sin(t)*cos(p));
define(y(r,t,p), r*sin(t)*sin(p));
define(z(r,t,p), r*cos(t));
depends(f, [x, y, z], x, [r, t, p], y, [r, t, p], z, [r, t, p]);
diff(f, r);
subst(diff(x(r, t, p), r), diff(x, r), %);
subst(diff(y(r, t, p), r), diff(y, r), %);
subst(diff(z(r, t, p), r), diff(z, r), %);
/* subst(x(r,t,p), x, %) などでは微分演算子の分母まで置換してしまってエラーとなる */
dfdr: 'diff(f, r) = %;
diff(f, t);
subst(diff(x(r, t, p), t), diff(x, t), %);
subst(diff(y(r, t, p), t), diff(y, t), %);
subst(diff(z(r, t, p), t), diff(z, t), %);
dfdt: 'diff(f, t) = %;
diff(f, p);
subst(diff(x(r, t, p), p), diff(x, p), %);
subst(diff(y(r, t, p), p), diff(y, p), %);
subst(diff(z(r, t, p), p), diff(z, p), %);
dfdp: 'diff(f, p) = %;

solve([dfdr, dfdt, dfdp], [diff(f, x), diff(f, y), diff(f, z)]);
trigsimp(%);
sol: ratexpand(%); 
/* 解けた。diff(f, x, 2) などはどうやって求める? これを (rat)subst するのは無理だった
   diff(f, r, 2) などから同様にしてゴリ押しすることは可能だろうが、大変だ。
   derivsubst も違う……

   f に対して置換するとよさそうだ。
   depends の影響を受けて、また x, y, z が出現してしまうのを防ぐため、まず g に変える */
depends(g, [r, t, p]);

dfdx: rhs(sol[1][1]);
subst(dfdx, f, dfdx);
subst(g, f, %);
ddfdxdx: ev(%, diff);

dfdy: rhs(sol[1][2]);
subst(dfdy, f, dfdy);
subst(g, f, %);
ddfdydy: ev(%, diff);

dfdz: rhs(sol[1][3]);
subst(dfdz, f, dfdz);
subst(g, f, %);
ddfdzdz: ev(%, diff);

/* 以下で Laplacian も出るが、vect パッケージのほうがはるかに簡単である */
trigsimp(ddfdxdx + ddfdydy + ddfdzdz), ratexpand;