LL用向前差分格式计算初边值问题.doc

偏微分方程

课程设计

学号：0683110

姓名：陆莉

指导老师：翟方曼
2010.01
一．题目
用向前差分格式计算如下热传导方程的初边值问题

已知其精确解为

二．理论
作为模型，考虑一维热传导方程：
…………（1.1）
其中是正常数，是给定的连续函数。
现在考虑第二类初边值问题的差分逼近：
初始条件：…………（1.2）
边值条件：，，………（1.3）
假设和在相应区域光滑，并且在满足相容条件，使上述问题有惟一充分光滑的解。
1.建立差分格式
（1）.区域剖分
取空间步长和时间步长，其中都是正整数。用两族平行直线和将矩形域分割成矩形网格，网格节点为。以表示网格内点集合，即位于开矩形的网点集合；表示所有位于闭矩形的网点集合；是网格界点集合。
其次，用表示定义在网点的函数，
（2）.微分方程的离散，建立相应差分格式
将方程在节点离散化，
，…………（1.4）
对充分光滑的解，由Taylor展式：
…………（1.5）
…………（1.6）
…………（1.7）
（1.5）移项得：
…………（1.8）
（1.6）（1.7）相加得：
…………（1.9）
将（1.8）（1.9）代入（1.4）得：
…………（1.10）
其中，

舍去，得到逼近（1.1）的向前格式差分方程：
，……（1.11）
其中，，
记

则由（1.4）
由（1.11）
显然，截断误差

（3）.边界条件

在本题中，，，，，
2.稳定性分析
用傅里叶方法对差分格式进行稳定性分析
以表示网比，将（1.11）改写成便于计算的形式：
（本题中）
以代入，得

消去，则知增长因子

由，得

即
只需

解得
所以向前差分格式的稳定性条件是

MATLAB程序

取，，则，满足稳定性条件
另取，，则，亦满足稳定性条件
另取，，则，亦满足稳定性条件

formatlong
a=2;
l=1;
T=1;
N=10;
M=400;
h=l/N;
to=T/M;
r=(a*to)/h^2;
forj=1:N+1
x(j)=(j-1)*h;
fork=1:M+1
t(k)=(k-1)*to;
u(j,k)=exp(x(j)+2*t(k));
end
end
u%求解精确解

forj=1:N+1
x(j)=(j-1)*h;
us(j,1)=exp(x(j));
end
fork=1:M+1
t(k)=(k-1)*to;
us(1,k)=exp(2*t(k));
us(N+1,k)=exp(1+2*t(k));
end
fork=2:M+1
forj=2:N
us(j,k)=r*us(j-1,k-1)+(1-2*r)*u(j,k-1)+r*us(j+1,k-1);
end
end
us%求解数值解

fork=1:M+1
forj=1:N+1
R(j,k)=abs(u(j,k)-us(j,k));
end
end
R%计算误差

Rmax=max(max(R))%求误差的最大值

图---精确解与数值解的比较：
x=0:0.1:1;
holdon
plot(x,u(:,M+1),'b');
plot(x,us(:,M+1),'y');
title('t=1，h=1/10，τ=1/400时精确解和数值解的比较')
text(0.05,21,'蓝：精确解');
text(0.05,20,'黄：数值解');
holdoff

图---取不同步长时的误差比较：
x=0:1/10:1;
y=0:1/20:1;
z=0:1/40:1;
holdon
plot(x,R(:,M+1),'b');
holdoff
M分别取10，20，40

四．表格及图表
部分结点处的精确解、数值解和误差绝对值（取，）
精确解数值解（3，1）（0.2，0）1.2214031.2214030（3，41）（0.2，0.1）1.4918251.4916251.994660e-004（3，81）（0.2，0.2）1.8221191.8218472.719067e-004（3，121）（0.2，0.3）2.2255412.2252053.359077e-004（3，161）（0.2，0.4）2.7182822.7190674.107891e-004（3，201）（0.2，0.5）3.3201173.3196155.018075e-004（3，241）（0.2，0.6）4.0552004.0545886.129183e-004（3，281）（0.2，0.7）4.9530324.9522847.486214e-004（3，321）（0.2，0.8）6.0496476.0487339.143683e-004（3，361）（0.2，0.9）7.3890567.38793911.168120e-0