偏微分方程程序附件Word文档格式.docx
《偏微分方程程序附件Word文档格式.docx》由会员分享，可在线阅读，更多相关《偏微分方程程序附件Word文档格式.docx（23页珍藏版）》请在冰豆网上搜索。
forn=1:
5/tao-1
form=1:
10/h-1
U（m,n）=u（m+1,n+1）;
U=U'
;
figure
（1）;
gridon;
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;
fornn=1:
formm=1:
uu（mm,nn）=（nn*tao）^2*sin（mm*h）;
UU=uu'
error=U-UU;
figure
（2）;
mesh（X,T,error）;
%mesh（X,T,UU）;
迎风格式误差曲面图'
figure（3）;
mesh（X,T,UU）;
古典解'
2、Lax-Wendroff格式
functionlax_w1（h,lamda）%Lax-Wendroff二阶精度显示格式
10/h
ifm==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差分格式
Lax-Wendroff格式差分曲面图'
error=U-UU;
Lax-Wendroff格式误差曲面图'
3、隐式中心格式
functionimplicit_1（h,lamda）%隐式格式，精度为：
时间一阶，空间二阶；
无条件稳定
time1=5/tao;
%时间为1时的时间层网格点
space1=10/h;
%空间为1时的时间层网格点
space1;
time1;
r=lamda/2;
%使用追赶法求每一时间层的节点值
k=1:
space1-1;
%初始化系数矩阵
a（k,k）=0;
a（1,1）=1;
a（1,2）=r;
fork=2:
space1-2%定义系数矩阵
a（k,k-1）=-r;
a（k,k）=1;
a（k,k+1）=r;
a（space1-1,space1-1）=1;
p=1:
q=1:
time1-1;
Z（p,q）=0;
%矩阵a*x（i）=Z（i）;
time1
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））;
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））;
%调用追赶法程序，程序在下面
forv=1:
u（v+1,n）=E（v）;
隐式格式差分曲面图'
error=UU-U;
隐式格式误差曲面图'
调用追赶法接三对角矩阵函数chase.m为：
functionX=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）;
fori=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）;
X（n）=y（n）/U（n）;
fori=n-1:
-1:
1
X（i）=（y（i）-c（i）*X（i+1））/U（i）;
二、计算下面双曲方程的近似解
1、二阶显示格式
functionexplicit_2（h,lamda）%波动方程显示格式，二阶精度，lamda<
1时稳定
tao=h*lamda;
1-tao;
1-h;
time1=1/tao+1;
space1=1/h+1;
time1%边界条件离散
u（space1,n）=sin（（n-1）*tao）;
form=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））;
forn=3:
1/tao
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）;
%波动方程显示格式
time1-2
space1-2
波动方程显示格式差分曲面图'
%计算古典解
space1-2;
time1-2;
uu（mm,nn）=sin（mm*h*nn*tao）;
波动方程显示格式误差曲面图'
波动方程古典解'
2、隐式二阶格式
functionimplicit_2（h,lamda）%波动方程隐式格式，二阶精度，无条件稳定
%离散化的网格点
v（i,j）=0;
%定义关于时间的一阶导数
fori=1:
v（i,1）=（i-1）*h;
%时间的一阶导数初始值离散化
forj=1:
w（i,j）=0;
%定义关于空间的一阶导数
end
w（i,1）=0;
%空间的一阶导数初始值离散化
%使用追赶法求每一时间层v的节点值,需要定义系数矩阵a和非其次项Z
r=-0.25*lamda^2;
m=1+0.5*lamda^2;
a（1,1）=m;
space1-3%定义系数矩阵
a（k,k-1）=r;
a（k,k）=m;
a（space1-2,space1-3）=r;
a（space1-2,space1-2）=m;
%定义非齐次项
%使用隐式格式计算
forjj=2:
time1%从第二时间层开始计算
v（1,jj）=0;
%v左边界值离散
v（space1,jj）=（u（space1,jj）-u（space1,jj-1））/tao;
%v右边界值离散
%计算非齐次项
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））;
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））;
jj-1））;
forkk=1:
v（kk+1,jj）=E（kk）;
%计算w（:
jj）和u（:
jj）
forii=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））;
forg=2:
u（g,jj）=u（g,jj-1）+tao*v（g,jj）;
i1=1:
uu（i1,jj）=0;
三、计算下面抛物方程的近似解
1、向前差分格式程序：
functiondiffusion_1（h,lamda,time）
%扩散方程向前差分格式，时间一阶精度，空间二阶精度，lamda<
=0,5稳定
%h为空间步长，lamda为网格比，t为演化时间
%h=0.1;
time=1;
lamda=0.5;
tao=h^2*lamda;
t=0:
time;
x=0:
1;
time1=time/tao+2;
wucha（1,i）=0;
space1-1%边界条件初始化
u（i,1）=sin（pi*（i-1）*h）;
u（m,n）=（1-2*lamda）*u（m,n-1）+lamda*u（m-1,n-1）+lamda*u（m+1,n-1）;
u=u'
mesh（X,T,u）;
向前差分格式曲面图'
uu=exp（-pi^2.*T）.*sin（pi.*X）;
error=uu-u;
向前差分格式误差曲面图'
%plot（error）;
mesh（X,T,uu）;
wucha=error（22,:
）
figure（4）
plot（wucha,'
ro-'
t=0.1s，x轴'
误差轴'
0.1秒时差分格式解与精确解的误差'
2、向后差分格式程序
functiondiffusion_2（h,lamda,time）
%扩散方程向后差分格式，时间一阶精度，空间二阶精度，无条件稳定
%使用追赶法求每一时间层u的节点值,需要定义系数矩阵a和非其次项Z
r=-lamda;
m=1+2*lamda;
kk=1:
a（kk,kk）=0;
time1%从第二时间层开始计算
Z（m,jj-1）=u（m+1,jj-1）;
%追赶法
u（kk+1,jj）=E（kk）;
向后差分格式曲面图'
向后差分格式误差曲面图'
figure（3）
0.1秒时向后差分格式解与精确解的误差'
3、Crank-Nicolson格式程序
functiondiffusion_3（h,lamda,time）
%扩散方程Crank-Nicolson差分格式，时间二阶精度，空间二阶精度，无条件稳定
r=-0.5*lamda;
r1=-r;
m=1+lamda;
m1=1-lamda;
%使用隐式格式计算
forjj=2:
Z（m,jj-1）=r1*u（m,jj-1）+m1*u（m+1,jj-1）+r1*u（m+2,jj-1）;
end
升级会员 