偏微分方程程序附件

程序

一、 计算下面双曲方程的近似解

解法：

1、 一阶迎风格式

程序：

function upwind1(h,lamda) %迎风格式 一阶精度 显示格式 x=h:h:10-h;

tao=lamda*h;

t=tao:tao:5-tao;

[X,T]=meshgrid(x,t);

i=1:10/h;

j=1:5/tao;

u(i,j)=0;

for n=2:5/tao

for m=2:10/h

u(m,n)=u(m,n-1)-lamda*(u(m,n-1)-u(m-1,n-1))+2*tao*(n-1)*tao*sin(m*h)+(n-1)^2*tao^3*cos(m *h); %迎风格式

end

end

for n=1:5/tao-1

for m=1:10/h-1

U(m,n)=u(m+1,n+1);

end

end

U=U';

figure(1);

grid on;

mesh(X,T,U);

xlabel('x 轴');

ylabel('t 轴');

zlabel('u 轴');

title('迎风格式差分曲面图');

%计算误差

⎪⎪⎩

⎪⎪⎨⎧==<<<<+=∂∂+∂∂== 0| 0|105,00,cos sin 2 0

02x t u u x t x t x t x u t u

jj=1:5/tao-1;

uu(ii,jj)=0;

for nn=1:5/tao-1

for mm=1:10/h-1

uu(mm,nn)=(nn*tao)^2*sin(mm*h);

end

end

UU=uu';

error= U-UU;

figure(2);

mesh(X,T,error);

%mesh(X,T,UU);

xlabel('x轴');

ylabel('t轴');

zlabel('u轴');

title('迎风格式误差曲面图');

figure(3);

mesh(X,T,UU);

xlabel('x轴');

ylabel('t轴');

zlabel('u轴');

title('古典解');

2、Lax-Wendroff格式

function lax_w1(h,lamda) %Lax-Wendroff 二阶精度显示格式

x=h:h:10-h;

tao=lamda*h;

t=tao:tao:5-tao;

[X,T]=meshgrid(x,t);

i=1:10/h;

j=1:5/tao;

u(i,j)=0;

for n=2:5/tao

for m=2:10/h

if m==10/h

u(m,n)=u(m,n-1)-lamda*(u(m,n-1)-u(m-1,n-1))+2*tao*(n-1)*tao*sin(m*h)...

+(n-1)^2*tao^3*cos(m*h); %用迎风格式处理边界问题

else

u(m,n)=u(m,n-1)-0.5*lamda*(u(m+1,n-1)-u(m-1,n-1))+2*tao*(n-1)*tao*sin(m*h)...

+tao*((n-1)*tao)^2*cos(m*h)+0.5*lamda^2*(u(m+1,n-1)-2*u(m,n-1)+u(m-1,n-1))... +tao^2*(n-1)*tao*(sin(m*h)-cos(m*h))+0.5*tao^2*(n-1)^2*tao^2*(sin(m*h)+cos(m*h)); %Lax-Wendroff差分格式

end

end

for n=1:5/tao-1

for m=1:10/h-1

U(m,n)=u(m+1,n+1);

end

end

U=U';

figure(1);

grid on;

mesh(X,T,U);

xlabel('x轴');

ylabel('t轴');

zlabel('u轴');

title('Lax-Wendroff格式差分曲面图');

%计算误差

ii=1:10/h-1;

jj=1:5/tao-1;

uu(ii,jj)=0;

for nn=1:5/tao-1

for mm=1:10/h-1

uu(mm,nn)=(nn*tao)^2*sin(mm*h);

end

end

UU=uu';

error=U-UU;

figure(2);

mesh(X,T,error);

%mesh(X,T,UU);

xlabel('x轴');

ylabel('t轴');

zlabel('u轴');

title('Lax-Wendroff格式误差曲面图');

3、隐式中心格式

function implicit_1(h,lamda) %隐式格式，精度为：时间一阶，空间二阶；无条件稳定x=h:h:10-h;

tao=lamda*h;

t=tao:tao:5-tao;

[X,T]=meshgrid(x,t);

time1=5/tao; %时间为1时的时间层网格点

space1=10/h; %空间为1时的时间层网格点