偏微分方程程序附件.doc

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

解法：
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('迎风格式差分曲面图');

%计算误差
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('迎风格式误差曲面图');

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
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时的时间层网格点
i=1:space1;
j=1:time1;
u(i,j)=0;
r=lamda/2;

%使用追赶法求每一时间层的节点值
k=1:space1-1; %初始化系数矩阵
a(k,k)=0;
a(1,1)=1;
a(1,2)=r;
for k=2:space1-2 %定义系数矩阵
a(k,k-1)=-r;
a(k,k)=1;
a(k,k+1)=r;
end
a(space1-1,space1-1)=1;

p=1:space1-1;
q=1:time1-1;
Z(p,q)=0; %矩阵a*x(i)=Z(i);

for n=2:time1
for m=2:space1-1
Z(m-1,n-1)=u(m,n-1)+tao*(2*(n-1)*tao*sin((m-1)*h)+...
n^2*tao^2*cos(m*h));
end
Z(space1-1,n-1)=u(space1,n-1)-lamda*(u(space1,n-1)-u(space1-1,n-1))...
+2*tao*(n-1)*tao*sin(space1*h)+...
(n-1)^2*tao^3*cos(space1*h); %边界离散，时间一阶，空间一阶
E=chase(a,Z(:,n-1)); %调用追赶法程序，程序在下面
for v=1:space1-1
u(v+1,n)=E(v);
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('隐式格式差分曲面图');

%计算误差
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=UU-U;
figure(2);
mesh(X,T,error);
%mesh(X,T,UU);
xlabel('x轴');
ylabel('t轴');
zlabel('u轴');
title('隐式格式误差曲面图');

调用追赶法接三对角矩阵函数chase.m为：

function X=chase(A,d)
%追赶法只用来解三对角方程组，其中A为方程的系数，d为右端项
%a为-1对角线上元素
%b为主对角线上的元素
%c为+1对角线上的元素
a=diag(A,-1);
b=diag(A);
c=diag(A,1);
n=length(b);
U(1)=b(1);
y(1)=d(1);
for i=2:n
L(i-1)=a(i-1)/U(i-1);
U(i)=b(i)-c(i-1)*L(i-1);
y(i)=d(i)-L(i-1)*y(i-1);
end
X(n)=y(n)/U(n);
for i=n-1:-1:1
X(i)=(y(i)-c(i)*X(i+1))/U(i);
end

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

解法：
1、 二阶显示格式

程序：

function explicit_2(h,lamda) %波动方程显示格式，二阶精度，lamda<1时稳定
tao=h*lamda;
t=tao:tao:1-tao;
x=h:h:1-h;
time1=1/tao+1; %时间为1时的时间层网格点
space1=1/h+1; %空间为1时的时间层网格点
[X,T]=meshgrid(x,t);
i=1:space1;
j=1:time1;
u(i,j)=0;

for n=2:time1 %边界条件离散
u(space1,n)=sin((n-1)*tao);
end

for m=2:1/h %初始条件离散 二阶精度
u(m,2)=u(m,1)+tao*h*(m-1)+0.5*lamda^2*(u(m+1,1)-2*u(m,1)+u(m-1,1));
end

for n=3:1/tao
for m=2:space1-1
u(m,n)=2*u(m,n-1)-u(m,n-2)+lamda^2*(u(m+1,n-1)-2*u(m,n-1)+u(m-1,n-1))+...
tao^2*(((n-1)*tao)^2-((m-1)*h)^2)*sin((n-1)*tao*(m-1)*h); %波动方程显示格式
end
end

for n=1:time1-2
for m=1:space1-2
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('波动方程显示格式差分曲面图');

%计算古典解
ii=1:space1-2;
jj=1:time1-2;
uu(ii,jj)=0;
for nn=1:time1-2
for mm=1:space1-2
uu(mm,nn)=sin(mm*h*nn*tao);
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、 隐式二阶格式

程序：

function implicit_2(h,lamda) %波动方程隐式格式，二阶精度，无条件稳定
tao=h*lamda;
t=tao:tao:1-tao;
x=h:h:1-h;
time1=1/tao+1; %时间为1时的时间层网格点
space1=1/h+1; %空间为1时的时间层网格点
[X,T]=meshgrid(x,t);
i=1:space1;
j=1:time1;
u(i,j)=0; %离散化的网格点
v(i,j)=0; %定义关于时间的一阶导数

for i=1:space1;
v(i,1)=(i-1)*h; %时间的一阶导数初始值离散化
end

for i=1:space1-1
for j=1:time1
w(i,j)=0; %定义关于空间的一阶导数
end
end
for i=1:space1-1
w(i,1)=0; %空间的一阶导数初始值离散化
end

for n=2:time1 %边界条件离散
u(space1,n)=sin((n-1)*tao);
end

%使用追赶法求每一时间层v的节点值,需要定义系数矩阵a和非其次项Z
r=-0.25*lamda^2;
m=1+0.5*lamda^2;
k=1:space1-2; %初始化系数矩阵
a(k,k)=0;
a(1,1)=m;
a(1,2)=r;
for k=2:space1-3 %定义系数矩阵
a(k,k-1)=r;
a(k,k)=m;
a(k,k+1)=r;
end
a(space1-2,space1-3)=r;
a(space1-2,space1-2)=m;
%定义非齐次项
p=1:space1-2;
q=1:time1-1;
Z(p,q)=0; %矩阵a*x(i)=Z(i);

%使用隐式格式计算
for jj=2:time1 %从第二时间层开始计算
v(1,jj)=0; %v左边界值离散
v(space1,jj)=(u(space1,jj)-u(space1,jj-1))/tao; %v右边界值离散
%计算非齐次项
for m=1:space1-3
Z(m,jj-1)=v(m+1,jj-1)+lamda*(w(m+1,jj-1)-w(m,jj-1))-...
r*(v(m+2,jj-1)-2*v(m+1,jj-1)+v(m,jj-1));
end
Z(space1-2,jj-1)=-r*v(space1,jj)+v(space1-1,jj-1)+lamda*...
(w(space1-1,jj-1)-w(space1-2,jj-1))-...
r*(v(space1,jj-1)-2*v(space1-1,jj-1)+v(space1-2,jj-1));
E=chase(a,Z(:,jj-1));
for kk=1:space1-2
v(kk+1,jj)=E(kk);
end
%计算w(:,jj)和u(:,jj)
for ii=1:space1-1
w(ii,jj)=0.5*lamda*(v(ii+1,jj)-v(ii,jj))+w(ii,jj-1)+...
0.5*lamda*(v(ii+1,jj-1)-v(ii,jj-1));
end
for g=2:space1-1
u(g,jj)=u(g,jj-1)+tao*v(g,jj);
end
end

for n=1:time1-2
for m=1:space1-2
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('波动方程显示格式差分曲面图');
%计算古典解
i1=1:space1-2;
jj=1:time1-2;
uu(i1,jj)=0;
for nn=1:time1-2
for mm=1:space1-2
uu(mm,nn)=sin(mm*h*nn*tao);
end
end
UU=uu';
error= U-UU;
figure(2);
mesh(X,T,error);
%mesh(X,T,UU);
xlabel('x轴');
ylabel('t轴');
zlabel('u轴');
title('波动方程显示格式误差曲面图');

三、 计算下面抛物方程的近似解

1、 向前差分格式程序：

function diffusion_1(h,lamda,time)
%扩散方程向前差分格式，时间一阶精度，空间二阶精度，lamda<=0,5稳定
%h为空间步长，lamda为网格比，t为演化时间
%h=0.1;time=1;lamda=0.5;

tao=h^2*lamda;
t=0:tao:time;
x=0:h:1;
[X,T]=meshgrid(x,t);
time1=time/tao+2; %时间为1时的时间层网格点
space1=1/h+1; %空间为1时的时间层网格点
i=1:space1;
j=1:time1;
u(i,j)=0; %离散化的网格点
wucha(1,i)=0;
for i=1:space1-1 %边界条件初始化
u(i,1)=sin(pi*(i-1)*h);
end
for n=2:time1
for m=2:space1-1
u(m,n)=(1-2*lamda)*u(m,n-1)+lamda*u(m-1,n-1)+lamda*u(m+1,n-1);
end
end
u=u';
figure(1);
grid on;
mesh(X,T,u);
xlabel('x轴');
ylabel('t轴');
zlabel('u轴');
title('向前差分格式曲面图');

%计算误差
ii=1:space1;
jj=1:time1;
uu=exp(-pi^2.*T).*sin(pi.*X);
error=uu-u;
figure(2);
mesh(X,T,error);
%mesh(X,T,UU);
xlabel('x轴');
ylabel('t轴');
zlabel('u轴');
title('向前差分格式误差曲面图');
%plot(error);
figure(3);
mesh(X,T,uu);
xlabel('x轴');
ylabel('t轴');
zlabel('u轴');
title('古典解');

wucha=error(22,:)
figure(4)
plot(wucha,'ro-');
xlabel('t=0.1s，x轴');
ylabel('误差轴');
title('0.1秒时差分格式解与精确解的误差');

2、 向后差分格式程序
function diffusion_2(h,lamda,time)
%扩散方程向后差分格式，时间一阶精度，空间二阶精度，无条件稳定
%h为空间步长，lamda为网格比，t为演化时间
%h=0.1;time=1;lamda=0.5;
tao=h^2*lamda;
t=0:tao:time;
x=0:h:1;
[X,T]=meshgrid(x,t);
time1=time/tao+2; %时间为1时的时间层网格点
space1=1/h+1; %空间为1时的时间层网格点
i=1:space1;
j=1:time1;
u(i,j)=0; %离散化的网格点
wucha(1,i)=0;
for i=1:space1-1 %边界条件初始化
u(i,1)=sin(pi*(i-1)*h);
end
%使用追赶法求每一时间层u的节点值,需要定义系数矩阵a和非其次项Z
r=-lamda;
m=1+2*lamda;
kk=1:space1-2; %初始化系数矩阵
a(kk,kk)=0;
a(1,1)=m;
a(1,2)=r;
for k=2:space1-3 %定义系数矩阵
a(k,k-1)=r;
a(k,k)=m;
a(k,k+1)=r;
end
a(space1-2,space1-3)=r;
a(space1-2,space1-2)=m;
%定义非齐次项
p=1:space1-2;
q=1:time1-1;
Z(p,q)=0; %矩阵a*x(i)=Z(i);
%使用隐式格式计算
for jj=2:time1 %从第二时间层开始计算
%计算非齐次项
for m=1:space1-2
Z(m,jj-1)=u(m+1,jj-1);
end
%追赶法
E=chase(a,Z(:,jj-1));
for kk=1:space1-2
u(kk+1,jj)=E(kk);
end
end
u=u';
figure(1);
grid on;
mesh(X,T,u);
xlabel('x轴');
ylabel('t轴');
zlabel('u轴');
title('向后差分格式曲面图');
%计算误差
uu=exp(-pi^2.*T).*sin(pi.*X);
error=uu-u;
figure(2);
mesh(X,T,error);
%mesh(X,T,UU);
xlabel('x轴');
ylabel('t轴');
zlabel('u轴');
title('向后差分格式误差曲面图');

wucha=error(22,:)
figure(3)
plot(wucha,'ro-');
xlabel('t=0.1s，x轴');
ylabel('误差轴');
title('0.1秒时向后差分格式解与精确解的误差');

3、 Crank-Nicolson格式程序
function diffusion_3(h,lamda,time)
%扩散方程Crank-Nicolson差分格式，时间二阶精度，空间二阶精度，无条件稳定
%h为空间步长，lamda为网格比，t为演化时间
%h=0.1;time=1;lamda=0.5;
tao=h^2*lamda;
t=0:tao:time;
x=0:h:1;
[X,T]=meshgrid(x,t);
time1=time/tao+2; %时间为1时的时间层网格点
space1=1/h+1; %空间为1时的时间层网格点
i=1:space1;
j=1:time1;
u(i,j)=0; %离散化的网格点
wucha(1,i)=0;
for i=1:space1-1 %边界条件初始化
u(i,1)=sin(pi*(i-1)*h);
end
%使用追赶法求每一时间层u的节点值,需要定义系数矩阵a和非其次项Z
r=-0.5*lamda;
r1=-r;
m=1+lamda;
m1=1-lamda;
kk=1:space1-2; %初始化系数矩阵
a(kk,kk)=0;
a(1,1)=m;
a(1,2)=r;
for k=2:space1-3 %定义系数矩阵
a(k,k-1)=r;
a(k,k)=m;
a(k,k+1)=r;
end
a(space1-2,space1-3)=r;
a(space1-2,space1-2)=m;
%定义非齐次项
p=1:space1-2;
q=1:time1-1;
Z(p,q)=0; %矩阵a*x(i)=Z(i);
%使用隐式格式计算
for jj=2:time1 %从第二时间层开始计算
%计算非齐次项
for m=1:space1-2
Z(m,jj-1)=r1*u(m,jj-1)+m1*u(m+1,jj-1)+r1*u(m+2,jj-1);
end
%追赶法
E=chase(a,Z(:,jj-1));
for kk=1:space1-2
u(kk+1,jj)=E(kk);
end
end
u=u';
figure(1);
grid on;
mesh(X,T,u);
xlabel('x轴');
ylabel('t轴');
zlabel('u轴');
title('Crank-Nicolson差分格式曲面图');
%计算误差
uu=exp(-pi^2.*T).*sin(pi.*X);
error=uu-u;
figure(2);
mesh(X,T,error);
%mesh(X,T,UU);
xlabel('x轴');
ylabel('t轴');
zlabel('u轴');
title('Crank-Nicolson差分格式误差曲面图');

wucha=error(22,:)
figure(3)
plot(wucha,'ro-');
xlabel('t=0.1s，x轴');
ylabel('误差轴');
title('0.1秒时Crank-Nicolson差分格式解与精确解的误差');

4、 Du Fort-Frankel格式程序
function diffusion_4(h,lamda,time)
%扩散方程Du Fort-Frankel差分格式，时间二阶精度，空间二阶精度，无条件稳定
%h为空间步长，lamda为网格比，t为演化时间
%h=0.1;time=1;lamda=0.5;
tao=h^2*lamda;
t=0:tao:time;
x=0:h:1;
[X,T]=meshgrid(x,t);
time1=time/tao+2; %时间为1时的时间层网格点
space1=1/h+1; %空间为1时的时间层网格点
i=1:space1;
j=1:time1;
u(i,j)=0; %离散化的网格点
wucha(1,i)=0;
for i=1:space1-1 %边界条件初始化
u(i,1)=sin(pi*(i-1)*h);
end
%使用追赶法求每一时间层u的节点值,需要定义系数矩阵a和非其次项Z
r=-0.5*lamda;
r1=-r;
m=1+lamda;
m1=1-lamda;
kk=1:space1-2; %初始化系数矩阵
a(kk,kk)=0;
a(1,1)=m;
a(1,2)=r;
for k=2:space1-3 %定义系数矩阵
a(k,k-1)=r;
a(k,k)=m;
a(k,k+1)=r;
end
a(space1-2,space1-3)=r;
a(space1-2,space1-2)=m;
%定义非齐次项
p=1:space1-2;
q=1;
Z(p,q)=0; %矩阵a*x(i)=Z(i);
jj=2; %使用Crank-Nicolson差分格式计算第二时间层
%计算非齐次项
for m=1:space1-2
Z(m,jj-1)=r1*u(m,jj-1)+m1*u(m+1,jj-1)+r1*u(m+2,jj-1);
end
%追赶法
E=chase(a,Z(:,jj-1));
for kk=1:space1-2
u(kk+1,jj)=E(kk);
end
%使用Du Fort-Frankel差分格式计算u值
d1=(1-2*lamda)/(1+2*lamda);d2=2*lamda/(1+2*lamda);
for j1=3:time1
for i1=2:space1-1
u(i1,j1)=d1*u(i1,j1-2)+d2*u(i1-1,j1-1)+d2*u(i1+1,j1-1);
end
end
u=u';
figure(1);
grid on;
mesh(X,T,u);
xlabel('x轴');
ylabel('t轴');
zlabel('u轴');
title('Crank-Nicolson差分格式曲面图');
%计算误差
uu=exp(-pi^2.*T).*sin(pi.*X);
error=uu-u;
figure(2);
mesh(X,T,error);
%mesh(X,T,UU);
xlabel('x轴');
ylabel('t轴');
zlabel('u轴');
title('Crank-Nicolson差分格式误差曲面图');

wucha=error(22,:)
figure(3)
plot(wucha,'ro-');
xlabel('t=0.1s，x轴');
ylabel('误差轴');
title('0.1秒时Crank-Nicolson差分格式解与精确解的误差');