XX学
科课程文
学 院: 理学院
专业年级: 09信息计算科学班
课 程: 偏微分方程数值解法
文题目:维热传导方程前Euler紧差分格式
指导教师:
2012年X月
学生姓名: 学 号:
分 工: 程序编写数值例子
学生姓名: 学 号:
分 工: 格式建立资料收集
学生姓名: 学 号:
分 工: 文档编辑资料整理
学生姓名: 学 号:
分 工: 公式编辑查找资料
学生姓名: 学 号:
分 工: 数分析查找资料
学生姓名: 学 号:
分 工: 数分析查找资料
目录
1引言 1
2物理背景1
3网格剖分 2
4.11前格式建立 2
412差分格式求解 4
413收敛性稳定性4
414 数值例子 7
421紧差分格式建立 10
422差分格式求解 12
423数值例子 13
总结17
参考文献 18
附录 19
1 引言
文考虑维非齐次热传导方程定解问题:
中正常数已知函数
目前常求解热传导方程差分格式前Euler差分格式Euler差分格式CrankNicolson格式Richardson格式.文出前Euler格式紧差分格式出截断误差数值例子.
2 物理背景
热传导物体部温度分布均匀热量物体温度较高点流温度较低点处函数表示物体时刻处温度假设关具二阶连续偏导数关具阶连续偏导数物体处热传导系数取正值设物体热容密度根热传导定律热量守恒定律公式
果物体均匀时均常数令式方程化
考虑物体热源热源密度函数热源热传导方程
中
3 网格剖分
取空间步长时间步长中正整数.两族行直线矩形域分割成矩形网格网格节点.记.表示网格点集合位开矩形网点集合表示位闭矩形网点集合网格界点集合
引进记号
分称穷范式(直范式)2范数(范数均范数)差商2范数(差商范式)范式
411前格式建立
定义网格函数
中
结点处考虑微分方程(311)
(32)
代入(32)
(331)
注意初边值条件(312)(313)
(332)
(333)
(331)~(333)中略量项
(34)
代差分格式
(351) (352)
(353)
称差分格式局部截断误差记
(361)
(362)
412 差分格式求解
记称步长差分格式(351)写
式表明第k+1层值第k层值显示表示出已知第k层值式直接第k+1层值时称(351)古典显格式古典显格式写成矩阵形式
413 收敛性稳定性
收敛性
设定解问题(311)~(313)解差分格式()()解时
中(361)定义
证明 记
(331)~(361)分(351)~(353)相减误差方程
时
证明完毕
稳定性
果应差分格式(351)~(353)时计算右端函数误差计算初值误差实际差分方程解
(38)
令
(351)~(353)(38)相减摄动方程组
(39)
时
(310)
式说明时误差
摄动方程组(39)差分方程(351)~(353)形式完全样述结果叙述
时差分格式(351)~(353)关初值右端项述意义稳定:设差分方程组
解
面考虑情况时必存时
设
(39)
验证解
易知
时
初始误差均会解较误差结
定理314 时差分格式(351)~(353)稳定
定理313定理314知步长时差分格式(351)~(353)稳定步长时差分格式(351)~(353)稳定种稳定性称条件稳定性实际计算时选取步长必须满足
414 数值例子
应前Euler 格式(351)(353)计算定解问题
述定解问题精确解
部分节点处数值解精确解误差绝值
数值解
精确解
|精确解数值解|
20
(0501)
18219e+000
18221e+000
23008e004
40
(0502)
20134e+000
20138e+000
34361e004
60
(0503)
22251e+000
22255e+000
41249e004
80
(0504)
24591e+000
24596e+000
46788e004
100
(0505)
27178e+000
27183e+000
52148e004
120
(0506)
30036e+000
30042e+000
57794e004
140
(0507)
33195e+000
33201e+000
63932e004
160
(0508)
36686e+000
36693e+000
70677e004
180
(0509)
40544e+000
40552e+000
78118e004
200
(05010)
44808e+000
44817e+000
86337e004
取长时数值解误差
110
1200
86337e004
*
120
1800
21748e004
39699e+000
130
13200
54366e005
40003e+000
140
112800
13591e005
40001e+000
误差曲面图()
误差曲面图()
误差曲线图
误差曲线
421紧差分格式建立
令
定义网格函数
Taylor展开
中式中标kk+1两等式相加2
利(342)
中
注意初值条件忽略项代差分格式
422 差分格式求解
差分格式写成
写出矩阵形式
中
时间层解三角线性方程组
423 数值例子
应紧差分格式(3471)(3473)计算定解问题
述定解问题精确解
误差
表中出空间步长缩原12时间步长缩原14时误差约缩原116
部分节点处数值解精确解误差绝值()
数值解
精确解
|精确解数值解|
10
(0501)
18221e+000
18221e+000
10836e006
20
(0502)
20138e+000
20138e+000
16247e006
30
(0503)
22255e+000
22255e+000
19547e006
40
(0504)
24596e+000
24596e+000
22195e006
50
(0505)
27183e+000
27183e+000
24750e006
60
(0506)
30042e+000
30042e+000
27435e006
70
(0507)
33201e+000
33201e+000
30351e006
80
(0508)
36693e+000
36693e+000
33554e006
90
(0509)
40552e+000
40552e+000
37087e006
100
(0510)
44817e+000
44817e+000
40990e006
取步长数值解误差
110
1100
40990e006
*
120
1400
25823e007
15873e+001
140
11600
16139e008
16000e+001
1800
14800
10090e009
15995e+001
1160
119200
65710e011
15355e+001
误差曲面图
误差曲面图
误差曲线图
误差曲线图
误差曲线图
总结:
文采差分格式紧差分格式求解抛物型方程差分格式采二层三点条件稳定显格式时误差着限增长稳定条件网出时误差紧差分格式采二层六点条件隐格式时间层均需解三角方程组
参 考 文 献
[1] 孙志忠.偏微分方程数值解法.北京:科学出版社2005
[2] 李荣华.偏微分方程数值解法.北京:高等教育出版社2005
[3]
附录A程序代码
流程图:
clear
ninput('n')空间剖分数
minput('m')时间剖分数
x(01n1)
t(01m1)
rn*nm网
l2
Adiag(ones(1n1)*(12*r))+diag(ones(1n2)*r1)+diag(ones(1n2)*r1)
精确解求解
for i1n+1
for j1m+1
u(ij)exp(x(i)+t(j))
end
end
初值条件求解
for i1n+1
u1(i1)exp(x(i))
end
for j1m+1
u1(1j)exp(t(j))
end
for j1m+1
u1(n+1j)exp(x(n+1)+t(j))
end
数值解求解
for j1m
u1(2nj+1)A*u1(2nj)+[r*u1(1j)zeros(1n3)r*u1(n+1j)]'
end
for i110
M(i)(u(6i*20+1))
N(i)(u1(6i*20+1))
end
abs(MN)'
误差图
uu1
[xt]meshgrid(xt)
surf(xtu'u1')
grid on
xlabel('x')
ylabel('t')
zlabel('|u(xt)u1(xt)|')
误差曲线
Eabs(u(2)u1(2))
plot(xE'o')
u(m+1)'u1(m+1)'
plot(xu(m+1)'u1(m+1)''o')
xlabel('x')
ylabel('|u(x1)u1(x1)|')
误差求解
for i2n
for j2m+1
K(i)max(u(ij)u1(ij))
end
end
Kmax(K)
附录B紧差分格式程序代码
clear
ninput('n')空间剖分
minput('m')时间剖分
H1n空间步长
T1m时间步长
x(0H1)
t(0T1)
rn*nm网
l2
Adiag(ones(1n1)*(56+r))+diag(ones(1n2)*(112r2)1)+diag(ones(1n2)*(112r2)1)k+1层系数矩阵
Bdiag(ones(1n1)*(56r))+diag(ones(1n2)*(112+r2)1)+diag(ones(1n2)*(112+r2)1)k层系数矩阵
精确解求解
for i1n+1
for j1m+1
u(ij)exp(x(i)+t(j))
end
end
初值条件求解
for i1n+1
u1(i1)exp(x(i))条件u(x0)exp(x)
end
for j1m+1
u1(1j)exp(t(j))条件u(0t)exp(t)
end
for j1m+1
u1(n+1j)exp(x(n+1)+t(j))条件u(1t)exp(1+t)
end
数值解计算
for j1m
u1(2nj+1)A\(B*u1(2nj)+[(112+r2)*u(1j)(112r2)*u(1j+1)zeros(1n3)(112+r2)*u(n+1j)(112r2)*u(n+1j+1)]')
end
for i1n
N(i)(u(6i*10+1))
M(i)(u1(6i*10+1))
end
abs(NM)'
误差曲面
uu1
[tx]meshgrid(tx)
surf(txu1u)
grid on
xlabel('x')
ylabel('t')
zlabel('|u(xt)u1(xt)|')
误差曲线
Eabs(u(m+1)'u1(m+1)')
plot(xE'o')
xlabel('x')
ylabel('|u(x1)u1(x1)|')
for i1n
M(i)(u(6i*10+1))
N(i)(u1(6i*20+1))
end
abs(MN)'
误差求解
for i2n
for j2m+1
K(i)max(abs(u(ij)u1(ij)))
end
end
Kmax(K)
中nan林业科技学
文档香网(httpswwwxiangdangnet)户传
《香当网》用户分享的内容,不代表《香当网》观点或立场,请自行判断内容的真实性和可靠性!
该内容是文档的文本内容,更好的格式请下载文档