《微分方程数值解》
课程论文
学生姓名1： 许慧卿 学 号： 20144329 学生姓名2： 向裕 学 号： 20144327 学生姓名3： 邱文林 学 号： 20144349 学生姓名4： 高俊 学 号： 20144305 学生姓名5： 赵禹恒 学 号： 20144359 学生姓名6： 刘志刚 学 号： 20144346 学 院： 理学院 专 业： 14级信息与计算科学 指导教师： 陈红斌
2017年6 月25日
《二阶常微分方程的差分格式》论文
一、《微分方程数值解》课程论文的格式 1、 引言：介绍研究问题的意义和现状 2、 格式：给出数值格式
3、 截断误差：给出数值格式的截断误差 4、 数值例子：按所给数值格式给出数值例子 5、 参考文献：论文所涉及的文献和教材
二、《微分方程数值解》课程论文的评分标准 1） 文献综述：10分；
2） 课题研究方案可行性：10分； 3） 数值格式：20分；
4） 数值格式的算法、流程图：10分； 5） 数值格式的程序：10分；
6） 论文撰写的条理性和完整性：10分； 7） 论文工作量的大小及课题的难度：10分； 8） 课程设计态度：10分； 9） 独立性和创新性：10分。
评阅人：
- 2 -
二阶常微分方程的差分格式
1 引言
考虑如下二阶常微分方程两点边值问题
d 2u (x ) -+u (x ) q (x ) =f (x ) , a
dx 2
u (a ) =α, u (b ) =β , (2) 的有限差分方法, 其中q (x ) , f (x ) 为[a , b ]上的连续函数, q (x ) ≥0; α, β为给定常数.
本文将对(1), (2)差分格式及紧差分格式进行求解, 并给出截断误差与数值例子.
2 差分格式及紧差分格式的建立 2.1 差分格式
将区间[a , b ]作M 等分, 记h =(b -a ) /M , x i =a +ih , 0≤i ≤M ,
Ωh ={x i |0≤i ≤M }, h 为网格步长, x i 为网格结点, Ωh 为网格, 称定义在Ωh 上的函数
为网格函数.
在网格结点x i 处考虑方程(1), 有
-u ''(x i ) +q (x i ) u (x i ) =f (x i ) , 1≤i ≤M -1, (3)
u (x 0) =α, u (x M ) =β.
任取一个结点x i , (i =1, 2, 3, … ,M -1) , 将u (x i +1), u (x i -1)在x i 处运用泰勒级数展开
u ''(x i ) h 2u '''(x i ) h 3u (4)(x i ) h 4
u (x i +1) =u (x i ) +u '(x i ) h ++++o (h 4) ,
2! 3! 4!
u '' (x i )(-h )u ''' (x i )(-h )u (
u (x i -1)=u (x i )+u ' (x i )(-h )+++
2! 3!
2
3
4)
(x i )(-h )
4!
4
+o (h 4).
利用上面两式得
u (4)(x i ) h 4
u (x i +1) +u (x i -1) =2u (x i ) +u ''(x i ) h ++o (h 4) .
12
2
整理得
u (x i +1) -2u (x i ) +u (x i -1) u (4)(x i ) h 22
''=u (x ) ++o (h ) . (4) i 2
h 12
由(4)式得
u ''(x u (x i +1) -2u (x i ) +u (x i -1) u (4)(x 2
i ) =h 2
-[i ) h 12
+o (h 2)]. (5) 将式(5)代入式(3)得
-
u (x i -1)-2u (x i )+u (x i +1)
h 2
+q (x i ) u (x i )=f (x i ) +R i .
(6)
其中R u (4)(x 2
i ) h i =-
12
+o (h 2) 为差分方程的截断误差. 用u i 代替u (x i ) , 并舍去截断误差得如下差分方程
-
u i -1-2u i +u i +1
h
2
+q (x i ) u i =f (x i ) . (7) 方程(7)结合边值条件, 得如下差分格式
-
u i -1-2u i +u i +1
h 2
+q (x i ) u i =f (x i ) , 1≤i ≤M -1. (8)
u 0=α, u M =β. (9)
2.2 紧差分格式
现在定义平均算子
⎧⎪1
12(u i -1+10u i +u i +1), 1≤i ≤M -1.
(Au ) ⎪⎪
i =⎨
⎪⎪
⎪⎩
u i , i =0, M .
在结点x i 处考虑方程(1), 有
-u ''(x i ) +q (x i ) u (x i ) =f (x i ) , 0≤i ≤M . (10)
用算子A 作用于(10)式, 整理可得
d 2u (x i ) d 2u (x i +1) ⎤11⎡d 2u (x i -1) 10
-⎢+10++q (x ) u (x ) +q (x i ) u (x i ) i -1i -1⎥22212⎣dx dx dx 12⎦12
+
根据引理
[1]
11
q (x i +1) u (x i +1) =[f (x i -1) +10f (x i ) +f (x i +1)].(14) 1212
6
, 当函数g (x ) ∈C [c -h , c +h ], 有
11h 4(6)
[g ''(c -h ) +10g ''(c ) +g ''(c +h ) ]=2[g (c +h ) -2g (c ) +g (c -h ) ]+g (ξi ) , 12h 240
ξi ∈(c -h , c +h ). (15)
结合(14), (15)两式, 整理有
-
11
u (x ) -2u (x ) +u (x ) +[q (x i -1) u (x i -1) +10q (x i ) u (x i ) +q (x i +1) u (x i +1)] []i -1i i +1h 212
1h 4(6)
=[f (x i -1) +10f (x i ) +f (x i +1) ]-u (ξi ) . (16) 12240
其中ξi ∈(x i -1, x i +1), 1≤i ≤M -1, u (x 0) =α, u (x M ) =β. 将(16)代入(10)中, 有
⎡u (x ) -2u (x i ) +u (x i +1) ⎤1
-⎢i -1+[q (x i -1) u (x i -1) +10q (x i ) u (x i ) +q (x i +1) u (x i +1)] 2⎥h ⎣⎦12
=
1
[f (x i -1) +10f (x i ) +f (x i +1) ]+R i , 12
h 4(6)
u (ξi ) 为差分方程的截断误差. 其中R i =-240
舍去截断误差, 并用u i 代替u (x i ) , 可得如下差分方程
-
u i -1-2u i +u i +11
+[q (x i -1) u i -1+10q (x i ) u i +q (x i +1) u i +1]
h 2121
=[f (x i -1) +10f (x i ) +f (x i +1) ],1≤i ≤M -1. (17) 12
方程(17)结合边值条件, 可得如下紧差分格式
u i -1-2u i +u i +11
+[q (x i -1) u i -1+10q (x i ) u i +q (x i +1) u i +1]
h 2121
=[f (x i -1) +10f (x i ) +f (x i +1) ], 1≤i ≤M -1, (18)
12-
u 0=α, u M =β. (19)
3 差分格式及紧差分格式的求解 3.1 差分格式的求解
将差分格式(8) , (9)中(8)整理为
-u i -1+[2+h 2q (x i )]u i -u i +1=f (x i ) h 2,
1≤i ≤M -1. (20)
可将式(20)表示为Au =F , 其中
⎛2+h 2q (x 1) ⎫-1 ⎪2
-12+h q (x ) -12 ⎪
⎪ A =
⎪2
-12+h q (x ) -1 ⎪M -2
-12+h 2q (x M -1) ⎪⎝⎭.
u =(u 1, u 2, , u M -2, u M -1), F =(f (x 1) +α,
即
T
f (x 2), , f (x M -2),
f (x M -1) +β),
T
⎛2+h 2q (x 1) ⎫⎛u 1⎫-1⎛ ⎪ ⎪ 2
u -12+h q (x ) -122 ⎪ ⎪
2
⎪ ⎪=h ⎪ ⎪ 2
u -12+h q (x ) -1 ⎪ M -2⎪M -2
2 ⎪ -12+h q (x M -1) ⎪⎝⎝⎭⎝u M -1⎭
f (x 1) ⎫⎛α⎫
⎪⎪f (x 2) ⎪ 0 ⎪⎪+ ⎪.
⎪ ⎪f (x M -2) ⎪ 0⎪
⎪f (x M -1) ⎪⎭⎝β⎭
该系数矩阵为三对角矩阵, 用追赶法进行求解较好. 3.2 紧差分格式的求解
将差分格式(18) , (19)中(18)整理为
⎡11⎤⎡11⎤⎡210⎤
-+q (x ) u ++q (x ) u +-+q (x ) j -1⎥j -1j ⎥j j +1⎥u j +1 ⎢⎢2⎢⎣h 12⎦⎣h 212⎦⎣h 212⎦
=
1101
f j -1+f j +f j -1, 1≤j ≤M -1. 121212
将上式整理为Au =F ,
其中
⎛ A =
⎝
210-11+q (x ) +q (x 2) 122h 12h 12-11210+q (x ) +q (x 3) 222h 12h 12
-11
+q (x M -2) 2
h 12
⎫⎪⎪⎪⎪ ，⎪⎪⎪210
+q (x ) M -1⎪h 212⎭(M -1)*(M -1)
u =[u 1, u 2, , u M -1],
T
⎛⎡11⎤⎫⎛10
⎢h 2-12q (x 0) ⎥α⎪ 12
⎦⎪ ⎣
⎪ 10 ⎪ 12F = ⎪+
⎪ 0 ⎪ ⎡1-1q (x ) ⎤β⎪
M ⎥ ⎢h 212⎪
⎣⎦⎭⎝⎝
即
1
1210
12
⎫⎛⎪ ⎪ ⎪ ⎪ ⎪1⎪
12⎪ 110⎪
⎪
1212⎭⎝f 1⎫⎛1⎫
f (x 0) ⎪⎪
f 2⎪12
⎪⎪ 0⎪⎪ ⎪ ⎪+ ⎪⎪ 0⎪⎪ ⎪f M -2⎪1
f (x M ) ⎪⎪⎭f M -1⎭⎝12
.
⎛
⎝
210-11+q (x ) +q (x 2) 122h 12h 12-11210+q (x ) +q (x 3) 2h 212h 212
-11
+q (x M -2) 2
h 12
⎫⎛u 1⎫⎪ ⎪⎪ u 2⎪⎪ ⎪⎪ ⎪ ⎪ ⎪⎪ ⎪⎪210 ⎪ +q (x ) ⎪M -1u M -1⎪2⎝⎭h 12⎭
⎫⎛f 1⎫⎛⎛⎡11⎤⎫⎛101
⎪ ⎪ ⎢h 2-12q (x 0) ⎥α⎪ 1212
⎦⎪ ⎪ f 2⎪ ⎣
⎪ ⎪ ⎪ 110 0
⎪ ⎪ ⎪ 1212
= ⎪ ⎪+ ⎪+ 1⎪
⎪ ⎪ 0
12⎪ ⎪ ⎪ ⎡1-1q (x ) ⎤β⎪ 110⎪ f M -2⎪
M ⎥ ⎢h 212⎪ ⎪ ⎪⎣⎦⎭⎝⎝1212⎭⎝f M -1⎭ ⎝
4 数值算例
算例1 应用差分格式(8), (9)计算如下两点边值问题
⎫
⎪⎪⎪0⎪
⎪. ⎪0⎪
⎪f (x M ) ⎪
⎪12⎭
f (x 0) 120
-u ''(x )+u (x )=e x (sin x -2cos x ), 0≤x ≤π, (21)
u (0)=0, u (π) =0. (22)
该定解问题的精确解为u (x )=e x sin x .
将区间[0,π]作M 等分, 并记h =π/M , x i =ih , 差分格式为
-u ''(x i )+u i =e x i (sin x i -2cos x i ), 1≤i ≤M -1 , u 0=0, u M =0.
下面将给出取不同步长时在结点处所得的最大模误差及阶
1) 最大模误差
E M =max |u i -u (x i ) |, (23)
i
2) 阶
rate =log 2
E M
. (25) E 2M
从下面表1可以看出, 步长缩小到原来的1/2时, 最大模误差缩小到原来的1/4. 下面图1给出了取不同步长时所得精确解与数值解的曲线.
表1 算例1取不同等分数目时最大模误差及其对应的阶
M
8 16 32 64 128
E M
1.532e-1 3.824e-2 9.555e-3 2.388e-3 5.971e-4
E M /E 2M
4.0004 4.0018 4.0005 4.0001 4.0000
rate
2.0023 2.0006 2.0017 2.0000 2.0000
图1 算例1取不同等分数目时数值解与精确解的对比曲线
由图1可知, 当等分数目M 逐渐增大时, 数值解与精确解的曲线拟合程度愈佳, 当M 增
大到一定程度时, 两者几乎重合。
算例2 应用差分格式(18),(19)计算算例1所给的两点边值问题(21),(22). 将区间[0,π]作M 等分, 并记h =π/M , x i =ih , 0≤x ≤π. 差分格式为
-
=
11u -2u +u +(u i -1+10u i +u i +1) i i +1)2(i -1h 12
1x i -1x i x i +1
⎡e sin x -2cos x +10e sin x -2cos x +e ()()(sin x i +1-2cos x i +1)⎤i -1i -1i i ⎦, 12⎣
1≤i ≤M -1 ,
u 0=0, u M =0.
下面表2给出了取不同步长时在结点处所得模误差及其阶. 可以看出, 当步长缩小到原来的1/2时, 最大模误差约缩小到原来的1/16. 下图2给出了算例2取不同步长时所得精确解与数值解的曲线.
表2 算例2取不同等分数目时最大模误差及其对应的阶
M
8 16 32 64 128
E M
1.874e-3 1.207e-4 7.583e-6 4.751e-7 2.970e-8
E M /E 2M
15.5312 15.9109 15.9609 15.9978 15.9947
rate
3.9571 3.9919 3.9965 3.9998 3.9995
图2 算例2取不同等分数目时的数值解与精确解的对比曲线
参考文献
[1] 孙志忠. 偏微分方程数值解法[M]. 第二版. 北京: 科学出版社, 2005.6 -7. [2] 李荣华. 偏微分方程数值解法[M]. 第二版. 北京: 高等教育出版社, 2005. [3] 刘卫国. MATLAB程序设计教程[M]. 第二版. 北京: 中国水利水电出版社, 2010.
◆ 程序流程图
总流程图及雅可比流程图如下：
◆ 程序代码
附录A
差分格式代码
建立函数文件Dirichlet.m, 如下：
function U=Dirichlet(M,h,qx,fx,u0,uM) %u0,uM为边值, M为区间[a0,b0]上的等分数目, h为步长, 其中qx, fx为(a,b)上的连续函数
A=zeros(M-1); %构造矩阵A 为(M-1)×(M-1) 的零矩阵
for i=1:M-1 %给矩阵A 主对角线上赋值
A(i,i)=2 + h*h*qx(i);
end
for j=1:M-2 %给矩阵A 对角线序号为1的位置赋值
A(j,j+1)=-1;
end
for k=2:M-1 %给矩阵A 对角线序号为-1的位置赋值
A(k,k-1)=-1;
end
B=zeros(M-1,1); %构造向量B 为(M-1)×1的零矩阵
for y=2:M-2
B(y,1)=h*h.*fx(y); %给向量B 从2到M-2行位置上赋值
end
B(1,1)=h*h.*fx(1)+u0; %给向量B 第一个元素赋值
B(M-1,1)=h*h.*fx(M-1)+uM; %给向量B 最后一个元素赋值
disp('数值解ux 为：'); % AU=B, 即 U=inv(A)*B
U=(A\B)'; %输出数值解的值
命令行窗口：
>> disp('输入定义区间[a,b]：')
输入定义区间[a,b]：
>> a=input('a=');
a=0
>> b=input('b=');
b=pi
>> M=input('输入等分数目M ：');
输入等分数目M ：8
>> h=(b-a)/M;
x=a+h:h:b-h;
qx=x./x;
ux=exp(x).*sin(x);
fx=exp(x).*(sin(x)-2.*cos(x));
>> u0=input('x=a处的函数值u(a)：')
x=a处的函数值u(a)：0
u0 =
>> uM=input('x=b处的函数值u(b)：')
x=b处的函数值u(b)：0
uM =
>> Dirichlet(M,h,qx,fx,u0,uM)
数值解ux 为：
Ans =
0.5303 1.4769 2.8905 4.6704 6.4287 7.3225 5.8940 >> disp('精确解ux 为：');
精确解ux 为：
>> ux
ux =
0.5667 1.5509 3.0009 4.8105 6.5819 7.4605 5.9796
附录B
紧差分格式代码
建立函数文件DirichletJ.m , 如下：
function U=DirichletJ(h,M,qx,fx,u0,uM,q0,qm,f0,fm) %u0,uM为边值, M 为区间[a0,b0]上的等分数目, h为步长, 其中qx, fx为(a,b)上的连续函数
A=zeros(M-1); %构造矩阵A 为(M-1)×(M-1)的零矩阵
for i=1:M-1 %给矩阵A 主对角线上赋值
A(i,i)=2./h^2 + 5*qx(i)./6;
end
for j=1:M-2 %给矩阵A 对角线序号为1的位置赋值
A(j,j+1)=-1./(h*h)+qx(j+1)./12;
end
for k=2:M-1 %给矩阵A 对角线序号为-1的位置赋值
A(k,k-1)=-1./(h*h)+qx(k-1)./12;
end
B=zeros(M-1,M-1); %构造向量B 为(M-1)×(M-1)的零矩阵
C=zeros(M-1,1); %构造向量C 为(M-1)×1的零矩阵
for c=1:M-1 %给矩阵B 主对角线上赋值
C(c,1)=fx(c);
end
for i1=1:M-1 %给矩阵B 主对角线上赋值
B(i1,i1)=5./6;
end
for j1=1:M-2 %给矩阵B 对角线序号为1的位置赋值
B(j1,j1+1)=1./12;
end
for k1=2:M-1 %给矩阵B 对角线序号为-1的位置赋值
B(k1,k1-1)=1./12;
end
D=B*C; %生成右边矩阵D
D(1,1)=D(1,1)+f0./12-u0./(h*h)-q0*u0./12; %给向量D 第一个元素赋值
D(M-1,1)=D(M-1,1)+fm./12-uM./(h*h)-qm*uM./12; %给向量D 最后一个元素赋值 disp('数值解ux 为：'); % AU=D, 即 U=inv(A)*D
U=(A\D)'; %输出数值解的值
命令行窗口：
>> disp('输入定义区间[a,b]：')
输入定义区间[a,b]：
>> a=input('a=');
a=0
>> b=input('b=');
b=pi
>> M=input('输入等分数目M ：');
输入等分数目M ：8
>> h=(b-a)/M;
x=a+h:h:b-h;
qx=x./x; %输入qx(x)表达式
q0=(a+1)./(a+1); %x=a处的函数值qx(a)
qm=(b+1)./(b+1); %x=b处的函数值qx(b)
ux=exp(x).*sin(x); %输入ux(x)表达式
fx=exp(x).*(sin(x)-2.*cos(x)); %输入fx(x)表达式
f0=exp(a).*(sin(a)-2.*cos(a));%x=a处的函数值fx(a)
fm=exp(b).*(sin(b)-2.*cos(b));%x=b处的函数值fx(b)
>> u0=input('x=a处的函数值u(a)：')
x=a处的函数值u(a)：0
u0 =
>> uM=input('x=b处的函数值u(b)：')
x=b处的函数值u(b)：0
uM =
>> DirichletJ(h,M,qx,fx,u0,uM,q0,qm,f0,fm);
数值解ux 为：
Ans=
0.5668 1.5509 3.0006 4.8097 6.5805 7.4586 5.9779 >> disp('精确解ux 为：');
精确解ux 为：
>> ux
ux=
0.5667 1.5509 3.0009 4.8105 6.5819 7.4605 5.9796