ストークス方程式のf

ストークス方程式の再現実験を行う時、いちいちfを計算するのは面倒だ。

maximaを使って自動化しよう。バッチファイルを作成します。

流れ：
u,p入力→f計算→ヘッダファイル出力

メモ帳で書いて拡張子を.macにしておきましょう.

kill(all)$
 
u:[sin(x+y+t),-sin(x+y+t)]$
pp:sin(x+y+t)$
 
diff(u[1],x)+diff(u[2],y);
p:pp-integrate(integrate(pp,y,0,1),x,0,1);
 
load(vect)$
f:diff(u,t)-ev(express(laplacian(u)),diff)+[diff(p,x),diff(p,y)];
 
load("cform.lisp")$
fp:openw("./func.h")$
printf(fp,"~a~%","double func_f1(double x,double y,double t){")$
printf(fp,"~t~a~a~a~%","return ",f[1],";")$
printf(fp,"~a~%","}")$
printf(fp,"~a~%","double func_f2(double x,double y,double t){")$
printf(fp,"~t~a~a~a~%","return ",f[2],";")$
printf(fp,"~a~%","}")$
printf(fp,"~a~%","double func_u1(double x,double y,double t){")$
printf(fp,"~t~a~a~a~%","return ",u[1],";")$
printf(fp,"~a~%","}")$
printf(fp,"~a~%","double func_u2(double x,double y,double t){")$
printf(fp,"~t~a~a~a~%","return ",u[2],";")$
printf(fp,"~a~%","}")$
printf(fp,"~a~%","double func_p(double x,double y,double t){")$
printf(fp,"~t~a~a~a~%","return ",p,";")$
printf(fp,"~a~%","}")$
close(fp)$
 

説明
kill(all)、メモリの初期化。奨励。
行末の"$"と";"、行末の記号。"$"は画面にアウトプットしない。";"はする。
u、想定する厳密解u。好きなように変える。
pp、想定する厳密解p。このあと積分平均0にするので、この時点で積分平均0になってなくてもよい。好きなように変える。
diff(u[1],x)+diff(u[2],y);、連続の式を満たすかチェック。結果を画面出力するための行末を";"にしている。
p、ppをもとに積分平均0にしたもの。積分値を引いてるだけ。
load(vect)、ラプラシアンやら勾配やら発散やらを計算するときにロードしておく。（なおこの状態で行列行列積を計算するとばぐるらしい）
f、ラプラシアンやら勾配やらを計算するところ。なおevとexpressがないと実際に計算してくれないので必要。
あとはcform.lispっていうc言語向けの出力を可能にするパッケージを入れて、ファイル出力。
なおprintfの書式については"maxima printf"でぐぐる。

maximaにパスを通してあればcygwinやらcmd上で

maxima -b func.mac
 

とすれば、以下のfunc.hが得られる.

double func_f1(double x,double y,double t){
 return 2*sin(y+x+t)+2*cos(y+x+t);
}
double func_f2(double x,double y,double t){
 return -2*sin(y+x+t);
}
double func_u1(double x,double y,double t){
 return sin(y+x+t);
}
double func_u2(double x,double y,double t){
 return -sin(y+x+t);
}
double func_p(double x,double y,double t){
 return sin(y+x+t)+sin(t+2)-2*sin(t+1)+sin(t);
}