再現実験で利用する

非定常の3 field Stokes方程式です

u,p,σを入力してそれらの導関数やf,gを出力します
pには項が加えられ自動的に積分平均0になります
3stokes.mac

kill(all);
 
u:[sin(x+y+t),-sin(x+y+t)];
pp:(x+y+t)^2;
sigma:[[cos(x+y+z)^2,x],[x,(x+y+z)^5]];
 
diff(u[1],x)+diff(u[2],y);
p:pp-integrate(integrate(pp,y,0,1),x,0,1);
integrate(integrate(p,y,0,1),x,0,1);
 
D:1/2*[[2*diff(u[1],x),diff(u[1],y)+diff(u[2],x)],[diff(u[2],x)+diff(u[1],y),2*diff(u[2],y)]];
f1:rho*diff(u[1],t)-2*etas*(diff(D[1][1],x)+diff(D[1][2],y))+diff(p,x)-(diff(sigma[1][1],x)+diff(sigma[1][2],y));
f2:rho*diff(u[2],t)-2*etas*(diff(D[2][1],x)+diff(D[2][2],y))+diff(p,y)-(diff(sigma[2][1],x)+diff(sigma[2][2],y));
g11:sigma[1][1]+lambda*diff(sigma[1][1],t)-2*etap*D[1][1];
g12:sigma[1][2]+lambda*diff(sigma[1][2],t)-2*etap*D[1][2];
g21:sigma[2][1]+lambda*diff(sigma[2][1],t)-2*etap*D[2][1];
g22:sigma[2][2]+lambda*diff(sigma[2][2],t)-2*etap*D[2][2];
 
nf:u[1]^2+u[2]^2+diff(u[1],x)^2+diff(u[2],x)^2+diff(u[1],y)^2+diff(u[2],y)^2+p^2+sigma[1][1]^2+2*sigma[1][2]^2+sigma[2][2]^2;
norm:integrate(integrate(nf,y,0,1),x,0,1);
 
u1x:diff(u[1],x);
u1y:diff(u[1],y);
u2x:diff(u[2],x);
u2y:diff(u[2],y);
 
load("cform.lisp");
with_stdout("f.txt",cform(u[1]),cform(u[2]),cform(u1x),cform(u1y),cform(u2x),cform(u2y),cform(p),cform(sigma[1][1]),cform(sigma[1][2]),cform(sigma[2][1]),cform(sigma[2][2]),cform(f1),cform(f2),cform(g11),cform(g12),cform(g21),cform(g22),cform(norm));
 
 

cform.lispは（どこかから入手して）C:\Program Files\Maxima-5.18.1\wxMaximaに置いてください。出力も同じフォルダに生成されます
wxMaximaを起動→ファイル→バッチ処理→3stokes.macを指定（数式が流れます）ダイバージェンスフリーか確認もできます
f.txt

sin(y+x+t);
-sin(y+x+t);
cos(y+x+t);
cos(y+x+t);
-cos(y+x+t);
-cos(y+x+t);
pow(y+x+t,2)-(6*pow(t,2)+12*t+7)/6.0;
pow(cos(z+y+x),2);
x;
x;
pow(z+y+x,5);
2*cos(z+y+x)*sin(z+y+x)+2*etas*sin(y+x+t)+rho*cos(y+x+t)+2*(y+x+t);
-5*pow(z+y+x,4)-2*etas*sin(y+x+t)-rho*cos(y+x+t)+2*(y+x+t)-1;
pow(cos(z+y+x),2)-2*etap*cos(y+x+t);
x;
x;
pow(z+y+x,5)+2*etap*cos(y+x+t);
(-495*cos(4*z+8)+990*cos(4*z+4)-7920*cos(2*z+4)+15840*cos(2*z+2)-495*cos(4*
    z)-7920*cos(2*z)+63360*pow(z,10)+633600*pow(z,9)+3326400*pow(z,8)[phone]*
    pow(z,7)[phone]*pow(z,6)[phone]*pow(z,5)[phone]*pow(z,4)[phone]*
    pow(z,3)[phone]*pow(z,2)[phone]*z-15840*cos(2*t+4)+31680*cos(2*
    t+2)-15840*cos(2*t)+42240*pow(t,2)+84480*t+2265904)/63360.0;
 
 

これをヘッダファイルにコピペします
func.h

double etas=1;
double etap=1;
double lambda=1;
double rho=1;
double alpha=1;
double uu1(double x,double y,double t){return sin(y+x+t);}
double uu2(double x,double y,double t){return -sin(y+x+t);}
double u1x(double x,double y,double t){return cos(y+x+t);}
double u1y(double x,double y,double t){return cos(y+x+t);}
double u2x(double x,double y,double t){return -cos(y+x+t);}
double u2y(double x,double y,double t){return -cos(y+x+t);}
double pp(double x,double y,double t){return pow(y+x+t,2)-(6*pow(t,2)+12*t+7)/6.0;}
double ss11(double x,double y,double t){return pow(cos(z+y+x),2);}
double ss12(double x,double y,double t){return x;}
double ss21(double x,double y,double t){return x;}
double ss22(double x,double y,double t){return pow(z+y+x,5);}
double ff1(double x,double y,double t){return 2*cos(z+y+x)*sin(z+y+x)+2*etas*sin(y+x+t)+rho*cos(y+x+t)+2*(y+x+t);}
double ff2(double x,double y,double t){return -5*pow(z+y+x,4)-2*etas*sin(y+x+t)-rho*cos(y+x+t)+2*(y+x+t)-1;}
double gg11(double x,double y,double t){return pow(cos(z+y+x),2)-2*etap*cos(y+x+t);}
double gg12(double x,double y,double t){return x;}
double gg21(double x,double y,double t){return x;}
double gg22(double x,double y,double t){return pow(z+y+x,5)+2*etap*cos(y+x+t);}
double maxima_norm(double t){return pow(z+y+x,5)+2*etap*cos(y+x+t);
(-495*cos(4*z+8)+990*cos(4*z+4)-7920*cos(2*z+4)+15840*cos(2*z+2)-495*cos(4*
    z)-7920*cos(2*z)+63360*pow(z,10)+633600*pow(z,9)+3326400*pow(z,8)[phone]*
    pow(z,7)[phone]*pow(z,6)[phone]*pow(z,5)[phone]*pow(z,4)[phone]*
    pow(z,3)[phone]*pow(z,2)[phone]*z-15840*cos(2*t+4)+31680*cos(2*
    t+2)-15840*cos(2*t)+42240*pow(t,2)+84480*t+2265904)/63360.0;}
 

メインプログラムにインクルードします

• 3stokes.c

#include"header.h"
#include"func.h"
#include"mat.h"
#include"mesh.h"
#include"solver.h"
#include"norm.h"
#include"from_maxima.h"
#include"output.h"
void apply_border_option(double **A,double *f,int *nl,int nn,int p1,double **nc,double t){
	int i,j;
	double tmp;
	for(i=0;i<nn;i++)
 
 
・・・
 
 
		free(u);
		free(AA);
	}
	fclose(fp);
	return 0;
}