3次元ラプラシアンの極座標表示をMaximaで求める
TwitterのTLで3次元ラプラシアンの極座標表示の手計算の話題が流れていました。手書きで解こうとは全く思えないのでフリーの数式処理ソフトのMaximaで解いてみます。数式処理ソフトを使う時点ですでに力づくなのですがなるべくさらっと見えるように書いてみます。
kill(all);
depends(f,[r,tht,phi]);
depends([r,tht,phi],[x,y,z]);
eqx:x=r*sin(tht)*cos(phi);
eqy:y=r*sin(tht)*sin(phi);
eqz:z=r*cos(tht);
eqxdx:diff(eqx,x);
eqxdy:diff(eqx,y);
eqxdz:diff(eqx,z);
eqydx:diff(eqy,x);
eqydy:diff(eqy,y);
eqydz:diff(eqy,z);
eqzdx:diff(eqz,x);
eqzdy:diff(eqz,y);
eqzdz:diff(eqz,z);
sol1:solve([eqxdx,eqxdy,eqxdz,eqydx,eqydy,eqydz,eqzdx,eqzdy,eqzdz],
[diff(r,x),diff(r,y),diff(r,z),diff(tht,x),diff(tht,y),diff(tht,z),
diff(phi,x),diff(phi,y),diff(phi,z)]);
sol1:trigsimp(sol1);
pol_lap:diff(f,x,2)+diff(f,y,2)+diff(f,z,2);
pol_lap:subst(sol1,pol_lap);
pol_lap:trigsimp(expand(pol_lap));
eqxd2x:diff(eqxdx,x);
eqxd2y:diff(eqxdy,y);
eqxd2z:diff(eqxdz,z);
eqyd2x:diff(eqydx,x);
eqyd2y:diff(eqydy,y);
eqyd2z:diff(eqydz,z);
eqzd2x:diff(eqzdx,x);
eqzd2y:diff(eqzdy,y);
eqzd2z:diff(eqzdz,z);
sol2:solve([eqxd2x,eqxd2y,eqxd2z,eqyd2x,eqyd2y,eqyd2z,eqzd2x,eqzd2y,eqzd2z],
[diff(r,x,2),diff(r,y,2),diff(r,z,2),diff(tht,x,2),diff(tht,y,2),diff(tht,z,2),
diff(phi,x,2),diff(phi,y,2),diff(phi,z,2)]);
sol2:trigsimp(sol2);
pol_lap:subst(sol2,pol_lap);
pol_lap:subst(sol1,pol_lap);
pol_lap:trigsimp(expand(pol_lap));
expand(%);
結果は
(%i34) expand(%);
2 2
d f d f
df ----- df -----
---- cos(tht) 2 2 -- 2 2
dtht dphi dr dtht d f
(%o34) ------------- + ------------ + ---- + ----- + ---
2 2 2 r 2 2
r sin(tht) r sin (tht) r dr
となって、ネット上の情報と比べると合っているようです。していることは
$${x=r sin(\theta) cos(\phi)}$$
$${y=r sin(\theta) sin(\phi)}$$
$${z=r cos(\theta)}$$
この座標変換の3式を$${x,y,z}$$で偏微分した式を偏微分項の連立方程式として解いて、直交座標のラプラシアンに代入していく。以上、終わり。
数式処理ソフトを使うと途中で式の項数が増えてきてもほぼ気にしなくていいので計算を簡単にするための工夫がいりません。当然でありますが変形の間違いや項のもれがないのもいいところです。ただし、変数の名前の入力を間違えていることはある(私の場合非常に多い)ので検算は手計算と同じように重要です。というか、単純な繰り返しはコード生成のコードを別言語で書いて自動生成するのが吉でしょう。私は(昔なつかしい)awkを使っています。今回は変数名を変えてコードを使いまわすわけでもないので直書きです。
この記事が気に入ったらサポートをしてみませんか?