NS方程由精确解求源项matlab代码

有时候在NS方程中我们需要构造有精确解的例子去测试格式的精度，这时候往往都是先给出精确解的表达式，然后带入原NS方程求出源项。如何求这个源项呢，对于可压流NS方程来说手动求导数会很麻烦。所以，这里我们用matlab符号计算来求源项。

clear all
syms x y u1 u2 u3 u4 pressure gamma txx txy tyx tyy mu Pr kTx kTy
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%各变量定义%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%x x方向坐标
%y y方向坐标
%u1 密度
%u2 x方向动量
%u3 y方向动量
%u4 总能
%pressure 压强
%gamma 气体常数
%txx txy tyx tyy 粘性应力张量
%mu 粘性系数
%Pr Prandtl数
%kTx $\kappa T_x$
%kTy $\kappa T_y$
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
u1 = sin(2*(x+y)) + 4;
u2 = sin(2*(x+y)) + 4;
u3 = sin(2*(x+y)) + 4;
u4 = (sin(2*(x+y)) + 4)^2;

pressure = (gamma - 1)*(u4 - 0.5*(u2^2+u3^2)/u1);
txx = 2.0/3.0 * mu * (2*diff(u2/u1,x)-diff(u3/u1,y));
txy = mu*(diff(u2/u1,y)+diff(u3/u1,x));
tyx = mu*(diff(u3/u1,x)+diff(u2/u1,y));
tyy = 2.0/3.0*mu*(2*diff(u3/u1,y)-diff(u2/u1,x));
kT = mu * gamma / Pr *(u4/u1-0.5*(u2^2+u3^2)/(u1*u1));
kTx = diff(kT,x);
kTy = diff(kT,y);

f1x = u2;
f1y = u3;
g1x = 0;
g1y = 0;

f2x = u2^2/u1 + pressure;
f2y = u2*u3/u1;
g2x = txx;
g2y = txy;

f3x = u2*u3/u1;
f3y = u3^2/u1 + pressure;
g3x = txy;
g3y = tyy;

f4x = u2/u1*(u4+pressure);
f4y = u3/u1*(u4+pressure);
g4x = u2/u1*txx + u3/u1*txy + kTx;
g4y = u2/u1*txy + u3/u1*tyy + kTy;

d1 = diff(f1x,x) + diff(f1y,y) - diff(g1x,x) - diff(g1y,y);
d2 = diff(f2x,x) + diff(f2y,y) - diff(g2x,x) - diff(g2y,y);
d3 = diff(f3x,x) + diff(f3y,y) - diff(g3x,x) - diff(g3y,y);
d4 = diff(f4x,x) + diff(f4y,y) - diff(g4x,x) - diff(g4y,y);

d1 = simplify(d1)
d2 = simplify(d2)
d3 = simplify(d3)
d4 = simplify(d4)