有限差分法求解偏微分方程Word文件下载.docx
- 文档编号:20647929
- 上传时间:2023-01-24
- 格式:DOCX
- 页数:35
- 大小:28.21KB
有限差分法求解偏微分方程Word文件下载.docx
《有限差分法求解偏微分方程Word文件下载.docx》由会员分享,可在线阅读,更多相关《有限差分法求解偏微分方程Word文件下载.docx(35页珍藏版)》请在冰豆网上搜索。
u(x,tk)
(_^)i,k(tki(汁)i,k
tk)o((tki
o
(2)
tk)2)
(2-4)
得到对时间的一阶偏导数
(U)=u(Xi,tki)u(Xi,tk)()i,k=
o()
(2-5)
由泰勒展开公式将u(x,t)对位置展开得
u(X,tk)u(Xi,tk)^U)i,k(XXi)
X
2
—)i,k(X
Xi)2
3
o((XXi))
当XXi禾口XXi
(2-6)
u(Xi,tk)
u(Xi,tk)
(_u)i,k(Xii
Xi)
2!
(
2)i,k(Xii
o((Xii
X))
u(N
i,tk)
(-^i'
kdii
2,(
2口
X)2
0((Ni
G3)
(2-7)
i时,代入式(2-6)得
因为XkiXkh,代入上式得
u(Xii,tk)
u(Xi,tk)
U)i,k
u(Xii,tk)
U(N,tk)
(」)i,k
)i,kh
2)i,kh
o(h3)
()i,k
o(h2)
(2-8)
得到对位置的二阶偏导数
将式(2-5)、(2-9)代入一般形式的抛物线型偏微分方程得
U(Xi,tki)u(x,tk)
u(xi,tk)2u(Xi,tk)u(Xii,tk)
h2
f(Xi,tk)o(h)
(2-10
)
为了方便我们可以将式
(2-10)
k1k
UiUi
写成
kkk
Ui12UiUi1
(2-11)
k
Ui
1k
声Ui1
2Uk
Ui1
(2-12)
最后得到古典显格式的差分格式为
k1
(12ra)ukr
(2-13)
其中r孑,古典显格式的差分格式的
截断误差是o(h2)o
2.1.2古典显格式稳定性分析
古典显格式(2-13)写成矩阵形式为
uh
12raI
raC
(2-14)
kzkkk
其中r—,Uh(U1,U2,……,Un
2,UN
1)
1
C0
M
0(N1)(N1)
上面的C矩阵的特征值是:
2cos(j
h)
H12ra
2rarac=
2ra2racos(jh)
2ra1cos(j
/-2j_±
4rasin
j1,2,……,N1
使(H)1,即
114rasin2」h1
c1
0ra-
结论:
当0ra
丄时,所以古典显格式是稳定的。
古典隐格式
2.2.1古典隐格式的推导
将ttk1代入式(2-3)得
U(Xj,tk1)U(Xj,tk)(~U)j,k(tk1tk)o((tk1
U(Xj,tk1)U(Xj,tk)(~U)j,k
(2-16)
(2-17)
。
()
(2-18)
U、U(Xj,tk)U(Xj,tk1)
;
)j,k=—
将式(2-9)、(2-18)
原方程得到
u(Xj,tk)u(Xj,tk1)
u(Xj1,tk)2u(Xj,tk)
u(Xj1,tk)
f(Xj,tk)o(h2)
Uj12UjUj1
(2-20)
(2-19
为了方便把(2-19)写成
kk1
UjUj
最后得到古典隐格式的差分格式为
kk
Uj12Uj
Uj
(2-21)
(12ra)u:
Uj1Uj1
(2-22)
其中r,古典隐格式的差分格式的截断误差是
h2)。
2.2.2古典隐格式稳定性分析
将古典隐格式(2-22)写成矩阵形式如下
12raIraCu;
1u;
f;
(r-j)
(2-23)
误差传播方程
12raIraCv;
1v;
(2-24)
A12raIraC,BI
所以误差方程的系数矩阵为
11
HA112raIraC
1
12ra2racosjh
使(H)1,显然
12ra2racos(jh)
12ra(1cos(jh))
14rasin2U
1恒成立。
对于r0,即任意网格比下,古典隐格式是绝对稳定的
Richardson格式
2.3.1Richardson格式的推导
将ttk1和ttk1,代入式(2-3)得
U2
(2-25)
u(冷tk1)u(Xi,tk)(t)i,k(tk1tk)o((tk1tk))
u(N,tk1)u(Xi,tk)(~u)i,k(tk1tk)o((tk1tk)2)
U(K,tki)U(N,tk)(「)i,ko()
U(X,tki)u(x,tk)(〒)i,ko()
(2-26)
由此得到可得
(2-27)
(<
)i,k
u(x,tk1)u(Xi,tk1)
将式(2-9)、(2-27)代入原方程得到下式
u(Xj,tkJu(x,tkJ
f(X,tQo(2h2)
u(Xii,tk)2u(x,tQu(Xii,tk)
(2-28)
为了方便可以把式(2-28)写成
k1k1
TTUi12UiUi1
(2-29)
(2-30)
最后得到Richardson显格式的差分格式为
k1ckckk
ui2rui12uiui1
2fik
(2-31)
其中r^2,古典显格式的差分格式的截断误差是o(2『)
2.3.2Richardson稳定性分析
将Richardson显格式(2-31)写成如下矩阵形式
k1ccIkk1c丄k
Uh2rC2Iuhuh2fh
误差传播方程矩阵形式
k1kk1
vh2rC2Ivhvh
VhVh
(2-32)
(2-33)
再将上面的方程组写成矩阵形式
Wh
vk12ra(C2I)I
vkI0
W
(2-34)
系数矩阵的特征值是
4racos(jh)4ra
j1
解得特征值为
28rasin2」h1
(2-35)
8rasin2-~~h
』64r2a2sin4~~h4
Y2
1,2
(2-36)
Max1
2=4rasin2
+Jl+16r2a2sin4~~-
(恒成立)
(2-37)
上式对任意的网比都恒成立,即Richardson格式是绝对不稳定的
4.Crank-Nicholson格式
3.4.1Crank-Nicholson格式的推导
将ttk1代入式(2-9)得
U(Xj,tk1)U(Xj,tk)(〒)j,k(tqtk)o((tk1tk))
2t22(2-40)
U(Xj,tk1)U(Xj,tk+1)([)j,k1(tkjtk+1)o((tk1tk+1))
即
U(Xj,tk1)U(Xj,tk)^^)j,k-o
(2)
2t2(2-41)
U(Xj,tk1)U(Xj,tk+1)(t)j,k1-o()
得到如下方程
u2u(Xj,tk;
)U(Xj,tk)2
(〒)j,k=0()
(2-42)
*)j,k1=
2u(Xj,tk2)u(Xj,tk1)
0
(2)
所以ttki处的一阶偏导数可以用下式表示:
2(寸)j,ki(
U-)j,k
12u(Xj,tk1)u(Xj,tk)2
u(Xj,tki)u(Xj,tki)
u(Xj,tki)u(Xj,tk)“2、
0()
(2-43)
将ttk,XXj1和X
Xj1代入式(2-6)可以得到式(2-9);
同理ttk1,XXj1和xXj1代入式(2-6)可以得到
(2-44)
uu(Xj1,tk1)2u(Xj,tk1)u(Xj1,tk1)亠2、
()j,k1~2o(h)
xh
所以ttk1,Xj处的二阶偏导数用式(2-6)、(2-44)表示:
_u)
2)j,k1
u(Xj1,tk1)2u(Xj,tk1)u(Xj1,tk1)u(Xj1,tk)2u(Xj,tk)u(Xj“tQ
2h2h2(2-45)
所以ttk1,Xj处的函数值可用下式表示:
2f(Xj,tk1)f(Xj,tk)
(2-46)
原方程变为:
(_^)j,k
()j,k1
2f:
1f:
(2-47)
将差分格式代入上式得:
U(Xj,tk1)U(Xj,tk)
u(Xj1,tk1)2u(Xj,tk1)u(Xj1,tk1)
u(Xji,tk)2u(Xj,tk)u(Xji,tk)
2f(Xj,tk1)f(Xj,tk)o(
h2)
(2-48)
为了方便写成:
k1kk
Uj1Uj12Uj
(2-49)
最后得到
Crank-Nicholson
格式的差分格式为
Uj1
1ru:
(2-50)
其中r,
格式的差分格式的
截断误差是o(
3.4.1Crank-Nicholson稳定性分析
Crank-Nicholson格式写成矩阵形式如下:
r
(1r)I-2
Uh=(1
Uh+
(2-51)
误差传播方程是:
(1
r)ITC
Vh
r)I
(2-52)
可以得到:
(1r)I-
r)1
;
C
(2-53)
(1r
)racos(j
(1r)racos(jh)
2jh
12rasin
2_
(2-54)
使(H)1即
12rasin2」-
21
12rasin2」-12rasin2」-
22
(2-55)
12rasin2」h
(2-56)
2jh2jh
12rasin12rasin
(2-57)
2rasin2rasin-
.2jh,.2jh
rasin1rasin-
(2-58)
上式恒成立。
Crank-Nicholson格式对任意网格比也是绝对稳定的
5.DuFortFrankie格式(Richardson格式的改进)将uik1(uik1uk)代入式(2-31)并化简得到DuFortFrankle
k1kk1k1kk1
ui2rui1uiuiui1ui
2fj
(2-59)
2ra)uk12ruik1u^(12ra)u:
(2-60)
可以证明DuFortFrankle
格式是绝对稳定的。
因为此格式是
Richardson格式
的改进格式,因此截断误差还是o(2h2)
总结
(1)古典显格式
古典显格式的差分格式为
uk1(12ra)Uikrukiukifik
截断误差:
o(h2)。
稳定性:
当0ra丄时,古典显格式稳定。
(2)古典隐格式
古典隐格式的差分格式为
(12ra)ukru:
1u:
1fjk
o(h2)o
任意网格比古典隐格式绝对稳定。
⑶Richardson显格式
Richardson显格式的差分格式为
k1ck/■*kkk1rk
Ui2rUi12uiUi1u2fi
o(2h2)稳定性:
任意网格比Richardson格式绝对不稳定
⑷Crank-Nicholson格式
Crank-Nicholson格式的差分格式为
/k1rk1k1」krkk1-k1-k
1rUjyUj1Uj1=1rUjyuj1uj1+-fjfj
o(2h2)。
Crank-Nicholson格式对任意网格比绝对稳定。
(5)DuFortFrankle格式
12ru^uik1(12ra)uk12fik
DuFortFrankle格式对任意网格比绝对稳定
三、MATLA实现五种差分格式对偏微分方程的求解及误差分析
精确数值解
上述偏微分方程的精确解是
2t
u(x,t)etsin(x)(0x1,t0)区域取值范围:
0x1,0t0.2。
用MATLAB寸精确解进行编程画出三维图像精确解程序如下:
closeall
clc
[x,t]=meshgrid(0:
:
1,0:
u=exp(-pi*pi*t).*sin(pi*x)
mesh(x,t,u);
surf(x,t,u);
xlabel('
x'
),ylabel('
t'
),zlabel('
u'
);
title('
精确数值解'
shadinginterp
图3-1精确数值解的三维图
(a)精确数值解X-Y平面(b)精确数值解X-Z平面
(c)精确数值解Y-Z平面
图3-2精确数值解的三个平面图
五种差分格式MATLABS序
3.2.1古典显格式
T=
X=
M=41
N=11
u=zeros(M,N);
%
ra=(T/(M-1))/((X/(N-1))A2);
%fprintf('
稳定性系数S=ra
disp(ra);
fori=2:
N-1
xx=(i-1)*(X/(N-1));
u(1,i)=sin(pi*xx);
end;
fork=1:
u(k,1)=0;
构造一个M行N列的矩阵用于存放时间网格比
为:
\n'
显示网格比
即t=0时刻赋值,边界条件
t和变量x
u(k,N)=0;
%x=0,x=1处的边界条件
M-1%矩阵是从y轴表示行k,x轴表示列i。
由行开始fori=2:
u(k+1,i)=((1-2*ra)*u(k,i)+ra*u(k,i+1)+ra*u(k,i-1));
%此处为古
典显格式
end
enddisp(u);
[x,t]=meshgrid(1:
N,1:
M);
%surf(x,t,u);
title('
古典显格式'
显示差分法求得的值
将区域划分成网格对每个点赋值再画图
此程序得到的是图3-3
图3-3古典显格式程序结果图(r=)
图3-4精确数值解、古典显格式程序结果的
Y-Z平面图(r=)
图3-5古典显格式在取不同网格比时的误差传播结果图
r=)红线表示精确解、
图3-6不同时间取值时精确解、与古典显格式的值对比图(网格比
蓝色线表示差分格式的解
图3-1、图3-3对比可以看出,精确解和古典显格式(网格比r=)的图形是致的。
图3-4精确数值解、古典显格式的Y-Z平面图结果可以看出古典显格式
的边界值和精确解一样。
图3-5是r分别取、、、时的误差传播图像,边界位置网
格数为5处的误差为得到的,可以看出r小于等于是稳定的;
但是r大于出现不
稳定现象。
很好的验证了古典显格式稳定性。
3.2.2
古典隐格式closeallclc
N=21
%网格比
fprintf('
稳定性系数S=ra为:
构造一个M行N列的矩阵用于存放时间t和变量x
%t=0
时刻的赋值,边界条件
u(1,i)=sin(pi*xx);
处的边界条件
隐格式的矩阵形式中的A矩阵赋值
%x=0,x=1
A=zeros(N-1,N-1);
A(1,1)=1+2*ra;
A(N-1,N-1)=1+2*ra;
A(1,2)=-ra;
A(N-1,N-2)=-ra;
form=2:
N-2
A(m,m-1)=-ra;
A(m,m)=1+2*ra;
A(m,m+1)=-ra;
%以下是——追赶法求u值
d=zeros(N-1,1);
%隐格式右边初始矩阵n=length(d);
U=zeros(n);
L=eye(n);
y=zeros(n,1);
x=zeros(n,1);
fori=1:
N-1d(i,1)=sin(pi*(i-1)*(N-1)));
end%隐格式右边矩阵赋值下循环对矩阵A进行LU分解
U(1,1)=A(1,1);
n
L(i,i-1)=A(i,i-1)/U(i-1,i-1);
U(i-1,i)=A(i-1,i);
%U的上次对角线即为A的上次对角线
U(i,i)=A(i,i)-L(i,i-1)*U(i-1,i);
M-1%外层大循环是求时间网格2,3,M的求解u
%以下是追赶法的求解过程
%追的过程即Ly=d的求解y
y(1,1)=d(1,1);
y(i,1)=d(i,1)-L(i,i-1)*y(i-1,1);
%赶的过程
x(n,1)=y(n,1)/U(n,n);
fori=n-1:
-1:
x(i,1)=(y(i,1)-U(i,i+1)*x(i+1,1))/U(i,i);
即Ux=y的求解x
%
u(k+1,i)=x(i,1)
d=x%
disp(u);
古典隐格式'
追赶法结束
更新右边矩阵每次外循环更换右边矩阵
显示差分法求得的值将区域划分成网格对每个点赋值再画图
此程序得到图是3-7
图3-7古典隐格式程序结果图
r=2)
图3-10不同时间取值时精确解、与古典隐格式的值对比图(网格比r=2)红线表示精确解、蓝色线表示差分格式的解
图3-7古典隐格式在r=2的图形与精确解是一致的。
图3-8精确数值解、古典隐格式的Y-Z平面图结果可以看出古典隐格式
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 有限 差分法 求解 微分方程