热传导方程几种差分格式的MATLAB数值解法比较(3)

end

DuFortFrankel格式：将（7）式与初始条件及第一类边界条件式联立并整理将格式改写为：

k+1k－1k

1+2r）U=（1－2r）U+2rU（j－1+jj

k 2rUkj+1+2τfj

（15）2，…，N－1，j=0，1，…，M－1） （k=1，

U0=?（j=0，1，…，M－1）jj k

U0=g1k，Uk2，…，N－1）L=g2k（k=1，

U

k+1 U 2 ? k+1 UM－2 k+1 UM－1

k+11

fk+rUk+Uk+1）

00

12 k

f2

+τ ?

fk M－2

rk+1

fk（UkM－1+M+UM） 2

（14）

由于此格式是3层格式，需要事先知道前面2个时间层上的解，第1层上的解通过对初始条件离散可得，第2层上的解使用前面的任意1种2层差分格式得到，再用此格式求解其余时间层上的解．程序核心代码如下：

%GeneratesecondrowusingimplicitschemeA=（1+2*r）*eye（M－2）；temp=－r*linspace（1，1，M－3）；A=A+diag（temp，1）+diag（temp，－1）；b=U（1，2：M－1）'；

b（1）=b（1）+r*U（1，1）；b（end）=b（end）+r*U（1，M）；

U（2，2：M－1）=A/b；%solveequationAU=b%GenerateremainingrowsofUforkt=3：N

forjx=2：M－1U（kt，jx）=（2*r*U（kt－1，jx－1）+2*r*U（kt－1，jx+1）+（1－2*r）*U（kt－2，jx））/（1+2*r）；

endend

此方程组的系数矩阵严格对角占优，差分方程组

解存在唯一．

程序如下：

functionU=HeatConductCrankNicolson（a，Init-Cond，LeftCond，RightCond，r，M，N，x，t）%InitializeparametersandUU=zeros（N，M）；F=zeros（N，1）；

%RightfunctionconditionsF=feval（f，（0：k：（N－1）*k）'）；%Generatefirstrow

U（1，：）=feval（InitCond，x）；%BoundaryconditionsU（：，1）=feval（LeftCond，t）；U（：，M）=feval（RightCond，t）；%GenerateremainingrowsofUforkt=2：N

A=（1+r）*eye（M－2）；

temp=－r/2*linspace（1，1，M－3）；A=A+diag（temp，1）+diag（temp，-1）；%B=（1-r）*eye（M-2）；

temp=r/2*linspace（1，1，M－3）；B=B+diag（temp，1）+diag（temp，－1）；b=B*U（kt－1，2：M－1）'；

b（1）=b（1）+r/2*（U（kt－1，1）+U（kt，1））；

b（end）=b（end）+r/2*（U（kt－1，M）+U（kt，M））；

U（kt，2：M－1）=A/b；%solveequationAU^（n+1）=b，b=BU^n

3数值实验

使用前面4种不同差分格式求解下面热传

导方程的初边值问题．

2

?u?u

=20＜t＜1，0＜x＜1?t?xu（x，0）=sin（πx）0≤x≤1u（0，t）=0，u（1，t）=00≤t≤1

2

t）=e－πtsin（πx）．此问题的真解为u（x，

{

定义误差e=u－U，

kkk

‖e‖2=‖u－U‖2=

12

（u（xj，tk）－Uk∑j）．Mj=1

百度搜索“77cn”或“免费范文网”即可找到本站免费阅读全部范文。收藏本站方便下次阅读，免费范文网，提供经典小说教育文库热传导方程几种差分格式的MATLAB数值解法比较(3)在线全文阅读。