抛物型方程的有限差分方法.docxVIP

抛物型方程的有限差分方法 一，求解问题 考虑一维非齐次热传导方程的定解问题 = /(x,r),0 xl,OtT……(1) ot ar w(x, 0) = °(/), 0 x / (2) w(0,z) = ?(/),班 1,/) = /?(r),0 tT (3) 其中 a为正长数, 为已知函数，0(O) = g(O),0(1) = 0(O), 式(2)为初值条件，(3)为边值条件。 二，网格剖分 取空间步长h = l/M和吋间步长T = T/ N,其中M、N都是整数。用两 族平行直线x = x{=加,(i = 0丄…,M)和t = tk = kr(i = 0丄…，N)将矩形 域G = {0x/;0?T}分割成矩形网格，网格结点为(兀厶)。以表示 网格内点集合，即位于开矩形G的网点集合；瓦表示所有位于闭矩形乙的网点 集合；几 ◎-乙是网格界点集合。 其次，用力表示定义在网点(兀?,切的函数，\iM-\^kN-\o 用适当的差商代替方程(1)中相应的偏微商。 三，差分格式 1, 向前差分 ?向前差分格式 +1 k k c k . k u.—比? u . — 2w. +仏］ =Q ——+Tf. t h2 1 Z =/(兀) 色0 =° =(p(x)M =UM }iM-\,0kN-\ 以r = az/h2为网比。将上式改写为便于计算的形式，则得以下向量形式 叮 =（1- 2r）u + 心二 + uM ） 4- r/（x., tk） liM-\,0kN-l 上式表示第k层的值显示表示出来。已知第k层的值则可以 直接得到第k+1的值仏严11 GV M｝。将上式表示成矩阵形式: （知1、 叭 ‘1一2厂 r （k、 终 ”(％/?) +曲、 知1 % r l-2r r k ■ ■ ■ —— ? ? ? ? ? ? ■ ■ ■ ? ■ ■ + ■ ■ ■ Jt+1 UM-2 r l-2r r k UM-2 k+l l“M-l ) J r l-2r? k \WA/-1 丿 /（兀如儿）+叱丿 ?向前Euler格式的截断误差 记 r du d2uyr ?+l k k c Ic , k u. 一 u. u.q - 2u. + u. t u h2 k 然后在幻进行泰勒展开从而得到 uk uk^ =u： +r^-\(x t 1 1 d才(兀 ? du - z 2 \ +厂乔?)+如) u uk+l -u. du d2u 2\ ?1（也）+厂丽1（也）+°（￡） T 显然得到截断误差 (w)= 4,v-[^i； TOC \o 1-5 \h \z 1 1 u . ? 「 = -r( )(--) +0(r+/2-) 12r 2 dr =0(T + h2) 其中（「）：是r 在矩形?_] XXj+i,tk t tk+l +某点的值。 3r 。厂 ?向前Euler格式的稳定性估计 向前Euler格式C = （l-2r）/ + rS ,其中，/是ml阶单位矩阵 0 ,C ,C 的特征值2； = l-2?* + 2rcos j7ch = l-4rsin2 5 = 为使|2；|14-Mt或一 5肉| = l-4rsir? 乎 51 + Me必须且只须 4 厂 sin?乎 W2 + M丿?二1,2,…,N-1 即4 r2,rl/20故向前差分格式当rl/2时稳定，当厂1/2时不稳定。 ?算例：应用向前Euler格式计算定解问题 du 32w =0,0 x l,0vf5 1, w(x,0) = e\Q x 1, w(0, r) = d, w(l, 0 = eu, 0 r 1 上述定解问题的精确解为讥兀r)=幺刊o 解：通过编程可求得不同步长比的数值解。表1给岀了取步长h=vio和 r = l/200 （步长比r=l/2）时计算得到的部分数值结果，数值解很好的逼近 精确解。表2给出了取步长h=l/10和r = 1 /100 （步长比r=l）时计算得到的 部分数值解。随着计算层数的增加，误差越来越人，数值结果无实用价值。 从图1JL取时间步长和空间步长都为10可以看到后面数值解和精确解的误差 非常去 前向差分格式不稳定。出表3给出了 “1/2时，取不同步长，数值解 的最大误差 Eg(h,t) = max | | 从表3可以看出当空间步长缩小到原来切川寸间步长缩小到原来的1A时，最 大误差约缩小到原来的1/4. 表1给出了取步长h=V10和 = 1 / 200吋计算得到的部分数值结果 k (x,t) 数值解 精确解 |精确解■数值解| 20 (0. 5, 0. 1) 1.822 1.8221 0.000166 40 (0. 5, 0. 2) 2.0134 2.0138 0.000303 60 (0. 5, 0. 3) 2. 2251 2.2