《计算物理学》课件第8章.ppt

1、第8章 有限差分方法 8.1 有关物理问题与数学物理方程8.2 有限差分原理8.3 矩形域中泊松方程的有限差分法8.4 差分方程的迭代解法 8.5 非矩形边界区域泊松方程的有限差分法 8.6 一维扩散方程的有限差分法 8.7 二维扩散方程的有限差分法8.8 一维波动方程的有限差分法 习题八 下篇下篇 计计 算算 物物 理理 学学 第8章 有限差分方法 8.1 有关物理问题与数学物理方程有关物理问题与数学物理方程 8.1.1 方程的导出方程的导出一个连续体，如气体、液体或固体，以及一个场，如电磁场、温度场等,其状态可用一时空函数u(x,y,z,t)来描写。在这个连续体或场中发生的物理过程遵循特定

2、的自然规律，其数学描述就是u随时空变化的方程，它们常以偏微分方程的形式出现。第8章 有限差分方法 1 描述稳定过程的泊松方程描述稳定过程的泊松方程(椭圆型方程椭圆型方程)例如在静电场中，其状态可用标势函数u(x,y,z)来描写，根据高斯定理有 (E)=(8.1)其中,为介质的介电常数，为自由电荷密度，。而电场强度E又可以表示为E=-u，将其代入上式可得-(u)=。对于均匀介质，为常数，于是有 (8.2)xyz u2第8章 有限差分方法 该方程为静电场中的泊松方程，其中表示又称为拉普拉斯方程。而在静磁场中，同静电场类似，有A=-j(8.3)该方程为静磁场中的泊松方程，其中A是磁矢势，是介质的磁导

3、率，j为传导电流密度。2222222xyz02 u2第8章 有限差分方法 2 描述输运过程的扩散方程描述输运过程的扩散方程(抛物型方程抛物型方程)以物质输运为例,设(x,y,z,t)为某一流体的密度函数，由于它的不均匀而发生扩散。按照扩散定律，物质的扩散流密度j正比于密度的梯度，两者反向，即有 j=-D，其中D为扩散系数。由质量守恒可导出连续性方程 (8.4)()Dt j第8章 有限差分方法 若D是均匀的，则有 (8.5)该方程为扩散方程。同样，在连续体中，由于温度T的不均匀而发生热传导，对于这种方程可用类似方程描写为 (8.6)其中K为温度传导率(假设均匀)，式(8.6)描述的方程即为热传导

4、方程。2DtTKtT2第8章 有限差分方法 3 描述振动传播过程的波动方程描述振动传播过程的波动方程(双曲型方程双曲型方程)对于交变电磁场，场强E、B 同场势A、u的关系为 (8.7)将上式代入Maxwell方程，可以得到场势满足的波动方程 (8.8)其中c=1/为波速。BAut AE22222222111ctuuct AAj第8章 有限差分方法 对于在连续介质中机械振动的传播过程，有如下类似的波动方程 (8.9)其中,u是介质的位移，r是杨氏模量，a是波速，f是外力。frtuau112222第8章 有限差分方法 8.1.2 方程的分类方程的分类以上我们按照数理方程所描述的物理过程将它们分成三

5、类，而在数学理论中，则按照这些方程的结构进行分类。二阶线性偏微分方程的一般形式可以写成(设自变量只有两个，即x和y)(8.10),(22222yuxuuyxFyuCyxuBxuA第8章 有限差分方法 1椭圆型方程椭圆型方程(B2-4AC0)如一维波动方程属于这一类型，方程(8.8)和(8.9)可以写成 (8.13)对比式(8.10)可知，B=0，A=1，C=-1。02222tuxu第8章 有限差分方法 8.1.3 边界条件和初始条件边界条件和初始条件从物理上讲，以上给出的数理方程仅适用于描写在一连续体或场的内部发生的物理过程。但仅靠这些方程还不足以完全确定物理过程的具体特征，因为物理过程的具体

6、特征还与连续体或场的初始状态和边界受到的外界影响有关。从数学上讲，一个偏微分方程会有无限多个解，偏微分方程加上边界条件和初始条件，才能构成一个定解问题。第8章 有限差分方法 1 第一类边界条件第一类边界条件若u代表方程中的未知函数，用表示方程适用区域D的边界。第一类边界条件为u|=u0(rb,t)(8.14)其中,u0(rb,t)是定义在上的已知函数，rb是相应边界点的位矢。在这种边界条件下边界上连续体或者场的状态是已知的。第8章 有限差分方法 2第二类边界条件第二类边界条件表达式为 (8.15)其中,n表示的外法线，q0(rb,t)是定义在上的已知函数。若u是电磁场的势，则代表场强；若u是密

7、度(或温度、或位移)，则代表流量(或热流，或应力)。当q0=0时称为第二类齐次边界条件。0(,)buqtnrnunu第8章 有限差分方法 3.第三类边界条件第三类边界条件表达式为 (8.16)000()(,)bua ubctnr第8章 有限差分方法 其中a0、b0 和c0是定义在上的已知函数。这一条件的物理情况较为复杂。对于热传导问题，它对应着连续体通过表面与外界发生辐射或对流等方式交换热量，表面热流正比于表面温度u与外界温度u0之差，即=k(u-u0)。对于给定的物理问题，其边界条件可能是很复杂的。例如在区域D的一部分边界1上有第一类条件，一部分边界2上有第二类条件，而其余部分则满足第三类条

8、件，或者是上述三式不能表示的另外的条件，需要具体分析边界的物理情况，然后予以确切的数学描述。nunu第8章 有限差分方法 4 初始条件初始条件与时间坐标t相联系，给出初始瞬间待求函数u在场域各处的值，即u|t=0=f1(r)(8.17)以及初始瞬间场域各处u对时间的变化率，即 (8.18)02()tuftr第8章 有限差分方法 通常仅含初始条件的定解问题称为初值问题(柯西问题)；没有初始条件而只有边界条件的定解问题称为边值问题；既有初始条件又有边界条件的定解问题，则称为混合问题(也称初边值问题)。对于由拉普拉斯方程构成的第一、第二和第三类边值问题，常分别称为狄利克雷、诺伊曼和洛平问题。第8章

9、有限差分方法 数学物理方程有许多是线性方程，对于这一类特定的物理问题往往有特定的解析方法，如分离变量法、积分变换法、格林函数法等等。解有时能用各种初等函数和超越的特殊函数来表达,但这些只限于比较典型的情况。更多的实际物理问题只能用非线性方程或方程组来描述，其求解方法更为复杂，只有少数问题有解析解,这时需要借助计算机来求偏微分方程的数值解。变分解法、有限元法、有限差分法、边界元等方法是数值解法中应用广、精度较高的几种常见解法，其共同特点是将求解偏微分方程的问题转化为求解线性代数方程组的问题。第8章 有限差分方法 8.2.1 差商公式差商公式将微分方程中的微商用差商代替(离散化)，是差分法求解偏微

10、分方程的基础。设u是坐标x的函数，取x=h等分坐标轴，结点坐标为xi，相应的有ui(i=1,2,)，则有泰勒展开(8.19)8.2 有限差分原理有限差分原理 iiiiiuhuhuhuu!3!2321第8章 有限差分方法(8.20)(8.21)(8.22)iiiiiuhuhuhuu!3!2321 iiiiiuhuhuhuu!3)2(!2)2(2322 iiiiiuhuhuhuu!3)2(!2)2(2322第8章 有限差分方法 1.误差为误差为O(h)的差商公式的差商公式由式(8.19)可得一阶向前差商公式 (8.23)由式(8.20)可得一阶向后差商公式 (8.24)hOhuuxuiii1dd

11、hOhuuxuiii1dd第8章 有限差分方法 式(8.21)-2式(8.19),可得二阶向前差商公式(8.25)式(8.22)-2式(8.20),可得二阶向后差商公式 (8.26)hOhuuuxuiiii212222dd hOhuuuxuiiii221222dd第8章 有限差分方法 2 误差为误差为O(h2)的差商公式的差商公式将式(8.25)代入式(8.19)，可得一阶向前差商公式 (8.27)将式(8.26)代入式(8.20),可得一阶向后差商公式 (8.28)212234ddhOhuuuxuiiii 212234ddhOhuuuxuiiii第8章 有限差分方法 式(8.19)-式(8.

12、20)/(2h),可得一阶中心差商公式 (8.29)式(8.19)式(8.20)/h2，可得二阶中心差商公式 (8.30)2112ddhOhuuxuiii 2211222ddhOhuuuxuiiii第8章 有限差分方法 8.2.2 差分格式的收敛性和稳定性差分格式的收敛性和稳定性用某种近似差商代替相应的微商，叫做一种差分格式。所谓差分格式的收敛性，是指当步长h0时，差分方程的解是否收敛于微分方程的解。显然只有那些具有收敛性的差分格式才是所需要的。第8章 有限差分方法 所谓差分格式的稳定性，是指误差u在运算过程中不会失控，即累计误差是否会无限增加。任何初值扰动对差分数值解的影响随时间推移不再增加

13、(强稳定)或在一段时间内有界(弱稳定)。有关差分法的数学理论对上述收敛性和稳定性都作了详细讨论。以下仅作结论性说明，不予以数学证明。第8章 有限差分方法 8.3.1 五点差分格式五点差分格式在二维情况下，泊松方程形式为 (8.31)若f(x,y)=0，则式(8.31)为拉普拉斯方程。8.3 矩形域中泊松方程的有限差分法矩形域中泊松方程的有限差分法),(22222yxfyuxuu第8章 有限差分方法 取x=y=h的正方形网格覆盖x-y平面(如图8.1所示)，结点坐标为(xi,yj)(i=1,2,N;j=1,2,M)。结点处的函数为u(xi,yj)=uij。在(i,j)点，利用中心差商公式(8.3

14、0)，则式(8.31)变为 (8.32)和22xu22yuijijjijijijiijfuuuuuhu411,1,1,122第8章 有限差分方法 式(8.32)中涉及到五个点，即中央点(i,j)及周围四点，如图8.1所示。这种差分格式的误差为O(h2)。差分法的数学理论已经证明，由上述五点格式得到的差分方程，在给定的边界条件下具有唯一的解，且当h0时趋于微分方程的解。第8章 有限差分方法 图8.1第8章 有限差分方法 8.3.2 矩形域的拉普拉斯方程矩形域的拉普拉斯方程 (8.33)【例例8.1】用有限差分法求解拉普拉斯方程，边界条件如图8.2所示。022222yuxuu第8章 有限差分方法

15、图8.2第8章 有限差分方法 若取h=5，如图8.2所示有三个内点，相应的u值记为u1、u2、u3。根据式(8.32)，可列出关于三个内点的差分方程组它的矩阵形式为0100404043232121uuuuuuu4 1 00 4 10 1 4321uuu1000 0 第8章 有限差分方法 解得。若取h=2.5，精度会提高，当然计算量会增大。有限差分法计算结果与本题精确解的比较如表8.1所示。786.26143.7786.1321uuu第8章 有限差分方法 表表8.1 有限差分法计算结果与精确解的比较有限差分法计算结果与精确解的比较第8章 有限差分方法 根据(8.32)式可以得到泊松方程的差分方程

16、组 (8.34)8.4 差分方程的迭代解法差分方程的迭代解法 ijijjijijijiijfuuuuuhu)4(11,1,1,122第8章 有限差分方法 该方程组中方程个数等于网格结点数，并且每个方程的左边最多有五项。故该方程组的系数矩阵具有大型、稀疏和带状的特点。为节省内存和机时，一般不采用消元法，而是采用迭代法，包括如下三种方法。(1)同步法。同步法。将式(8.34)化为 (8.35)12,11,11,1 ()4kkiji jiji jijijuuuuuh f第8章 有限差分方法 上式中右边的u用第k步的值，则左边得到第k+1步的值。这种迭代方法称为同步法，它需要两套内存，一套存第k步的值

17、，一套存第k+1步的值，收敛较慢。(2)异步法。由于u的计算是按照i、j 由小到大进行的，迭代到某一步在计算新的uij时，新的ui，j-1和ui-1，j已被算出，所以式(8.35)可以改写为 (8.36)这种迭代方法称为异步法，它只需一套内存，收敛较快。1112,11,11,1()4kkkkkiji jiji jijijuuuuuh f第8章 有限差分方法(3)逐次超松弛迭代法(SOR法)。在异步法的基础上，加权平均，得到SOR法。引入超松弛因子，将式(8.36)改写为 (8.37)1112,11,11,()(1)4kkkkkkiji jiji jijijijuuuuuh fu第8章 有限差分

18、方法 当=1时，式(8.37)就是式(8.36)。当在(1,2)内取一适当值时，可获得较快的收敛速度。如在例8.1中，当h=2.5并要求|u(k+1)-u(k)|10-3时，取不同的值，所需的迭代次数k也不同，如表8.2所示。第8章 有限差分方法 表表8.2 第8章 有限差分方法 从表8.2可以看出，值的选取对迭代次数k影响较大。显然，并非越大或越小时，迭代次数k越小。那么如何选取最佳的值使得迭代次数k越小呢？通常采用如下方法:(1)对于矩形区域，常数边界条件下超松弛因子常取为 (经验公式)1,11222/122mnmn第8章 有限差分方法 其中，n+1、m+1分别为x轴和y轴上的等分结点数。

19、而对于正方形区域，取为其中n+1为每边上的等分结点数。nsin12第8章 有限差分方法(2)在(1,2)内，依次取几个值，用少数几步迭代，找出最佳值。SOR法迭代过程中除了超松弛因子的选取外，关于初值的选取也很重要。事实上，关于上述各种迭代法的收敛性问题，可以证明，当h0或k时，差分方程的解将趋近于微分方程的解。收敛性与初值无关，但收敛速度与初值有关。关于初值的选取问题，若 f(x,y)=0，可取u|平均值为初始值，或任取初值。若f(x,y)0，可取u=0为初值。第8章 有限差分方法 另外，关于收敛的判别方法，可事先指定出误差的范围，当时，迭代终止。采用有限差分方法求解偏微分方程时，还需注意对

20、边界条件的处理。对于规则矩形边界而言，研究区域的网格结点落在边界上，第一类边界条件显然无需再做处理，而对于第二、三类边界条件，可采用以下差分格式。根据以下第三类边界条件 (8.38)()1(maxkijkijuu),(yxcnubau第8章 有限差分方法 为找到误差为O(h2)的边界条件的差分式，在图8.3所示的xd和yd处用式(8.27)中的一阶向前差商代替该处的，而在xu和yu处用式(8.28)中的一阶向后差商代替该处的，解出边界值，其结果为 (8.39)nu /nu /bhabubuyxhcubhabubuyxhcubhabubuyxhcubhabubuyxhcuimimmiimjnjn

21、jnnjiiiijjjj324),(2324),(2324),(2324),(2212132113211第8章 有限差分方法 图8.3第8章 有限差分方法 一般情况下，a和b是x、y的函数，上式应作相应变化。显然，当a=0时，上式即为第二类边界条件的差分格式。通常边界是不规则的，边界条件处理方法如下：对于第一类边界条件u|=u0，为边界。若结点在边界上，直接代入。若结点不在边界上，采用不对称网格方法，如图8.4所示，边界与网格线交点为D和B，对于泊松方程,则有 (8.40)2)1()1(2)(21)(2111122PABPABuuuhahhauahuahuxu第8章 有限差分方法 (8.41)

22、其中，。此时泊松方程在边界上的差分格式为 (8.42)2)1()1(2222PCDuuuhyu,1haha2PPDCBAfhuuuuu221)11()1(1)1(1第8章 有限差分方法 对于第二类边界条件，若结点在边界上(如图8.5所示)，则有 (8.43)0unucoscosuuuuuunxyxy ni nj n第8章 有限差分方法 其中，。其差分格式为 (8.44),huuxuPBhuuyuPD0cos)(cos)(huuuuuPDPB第8章 有限差分方法 图8.4 第8章 有限差分方法 图8.5 第8章 有限差分方法 若结点不在边界上，过P点向边界作垂线，交边界于P。用P的法向单位矢n做

23、为P的法向矢，实质上是把P近似为P，差分格式处理同式(8.44)。对于第三类边界条件的处理方法是前两类边界条件处理方法的组合。0)(unubau第8章 有限差分方法【例例8.2】试采用有限差分法编程求解拉普拉斯方程在网格结点处的值。222204030,(04,03)(3),0sin,04xxyyuuxyxyuy yuuxu第8章 有限差分方法 图8.6第8章 有限差分方法 如图8.6所示，取h=1，=1.25，迭代终止误差为=10-4。计算程序:program main integer i,j real x(5),y(4),u(5,4),t,p,q real,parameter:w=1.25,

24、eps=1e-4,h=1,pi=3.1415926 !w，eps分别表示和第8章 有限差分方法 open(1,file=example.dat)do i=1,5 !输入边界条件x(i)=(i-1)*hu(i,1)=sin(pi*x(i)/4)u(i,4)=0.0 enddo do j=1,4y(j)=(j-1)*hu(1,j)=y(j)*(y(j)-3)u(5,j)=0.0第8章 有限差分方法 enddo do i=2,4 !SOR法赋初值 do j=2,3u(i,j)=0.0enddo enddo10 p=0.0 do i=2,4do j=2,3t=u(i,j)第8章 有限差分方法 u(i,

25、j)=w*(u(i,j-1)+u(i-1,j)+u(i,j+1)+u(i+1,j)/4+(1-w)*u(i,j)q=abs(u(i,j)-t)if(q.gt.p)p=q!p为每次迭代的最大误差 enddo enddo if(p.gt.eps)goto 10 第8章 有限差分方法 do i=1,5!输出结果 do j=1,4write(*,*)i,j,u(i,j)write(1,*)i,j,u(i,j)enddoenddoend第8章 有限差分方法 计算结果:i j u(i,j)1 1 0.000000e+00 1 2 -2.000000 1 3 -2.000000 1 4 0.000000e+

26、00 2 1 7.071068e-01 2 2 -4.403373e-01 2 3 -6.375487e-01 2 4 0.000000e+00 3 1 1.000000第8章 有限差分方法 3 2 1.690345e-01 3 3 -1.098504e-01 3 4 0.000000e+00 4 1 7.071069e-01 4 2 2.263135e-01 4 3 2.911735e-02 4 4 0.000000e+00 5 1 0.000000e+00 5 2 0.000000e+00 5 3 0.000000e+00 5 4 0.000000e+00第8章 有限差分方法 8.5.1

27、圆形域中泊松方程的有限差分解圆形域中泊松方程的有限差分解1泊松方程的差分格式泊松方程的差分格式这里仍以二维泊松方程为例，该方程表示为2u=f(r,)(0rM),A(N,3),B(M,3)输入D,a0,b0,tmax,h,(a,b)1,2,3,4，P=2D/h2,P1=1/P+1,P2=1/P-1，N=1+a0/h，M=1+b0/h，kmax=tmax/。第8章 有限差分方法 原始的和消元后的系数矩阵A和B(若与t有关，应置于k的循环之外)call coef(N,a1,b1,a2,b2,h,P1,A)call coef(M,a3,b3,a4,b4,h,P2,B)do 8 k=1,kmax,2沿y

28、方向算。t=k*do i=2,N-1 do j=1,M第8章 有限差分方法 Rj=ui-1,j+2*P2*uij+ui+1,j R1=R1*b3+2*h*c3(xi,t)Rm=Rm*b4+2*h*c4(xi,t)call solve(M,B,R)do j=1,M Vij=Rj边界值和上的u值为 do j=1,M V1j=2*h*c1(yj,t)+4*b1*V2j-b1*V3j/(2*h*a1+3*b1)第8章 有限差分方法 Vnj=2*h*c2(yj,t)+4*b2*Vn-1,j-b2*Vn-2,j/(2*h*a2+3*b2)沿x方向算。t=t+do j=2,M-1 do i=1,N Ri=V

29、ij-1+2*P2*Vij+Vij+1 R1=R1*b1+2*h*c1(yj,t)Rn=Rn*b2+2*h*c2(yj,t)call solve(N,A,R)do i=1,N第8章 有限差分方法 uij=Ri边界和上的u值为 do i=1,N ui1=2*h*c3(xi,t)+4*b3*ui2-b3*ui3/(2ha3+3b3)uim=2*h*c4(xi,t)+4*b4*ui,m-1-b4*ui,m-2/(2*h*a4+3*b4)8 输出t,u 需要说明的是，子程序coef和solve与8.6节的相同。第8章 有限差分方法 8.7.4 二维显式格式二维显式格式对于式(8.76)，在(i,j,k

30、)点，用向前差商，用中心差商代替，则可以得到 (8.83)tu和22xu22yukijijijjiijjikijkijuuuuuuhDuu)2()2()(1111121第8章 有限差分方法 引入P=D/h2，可把式(8.83)化为 (8.84)利用式(8.84)可由k时的u直接求k+1时的u，不必解联立方程，故称显式格式。误差为O()+O(h2)。显式格式虽然简单，但它的稳定性是有条件的，要求(8.85)11111 ()(1 4)()kkijijijijijijuP uP uuuu14P 第8章 有限差分方法 此外，在设计程序时需用两套u的数组，一套存k时刻的，另一套存k+1时刻的。对于边界条

31、件式(8.80)，将用误差O(h)的向前(、边界)或向后(、边界)差商代替，就得到边界条件的差分式。nu第8章 有限差分方法 计算程序流程为 dimension x(N),y(M),u(N,M),V(N)输入D,a0,b0,tmax,h,P,(a,b,c)1,2,3,4，=h2P/D,N=1+a0/h,M=1+b0/h,kmax=tmax/。结点坐标与初值i=1,N,xi=(i-1)h；j=1,M,yj=(j-1)h，uij=u0(xi,yj)第8章 有限差分方法 时间演化 do 40 k=1,kmax t=k*do 10 i=2,N-1 do 10 j=2,M-1 10 Vij=(1-P)*

32、uij+P*(ui,j-1+ui,j+1+ui-1,j+ui+1,j)do 20 i=2,N-1 do 20 j=2,M-120 uij=Vij 第8章 有限差分方法 边界值 do 25 j=1,M u1j=(b1*u2j+h*c1)/(h*a1+b1)边界25 unj=(b2*un-1,j+h*c2)/(h*a2+b2)边界 do 30 i=1,N ui1=(b3*ui2+h*c3)/(h*a3+b3)边界30 uim=(b4*ui,m-1+h*c4)/(h*a4+b4)边界40 输出t,u 第8章 有限差分方法 8.8.1 显式差分格式显式差分格式考虑如下一维波动方程的定解问题 (8.86

33、)8.8 一维波动方程的有限差分法一维波动方程的有限差分法 2220max2212001112220(,),0,0(,)(,)(),(),0,ttuucf x txatttxu x tu x tR xR xtuaubcxnua ubcxan 第8章 有限差分方法 取x=h，t=进行离散化，节点坐标为 (8.87)max,2 ,1 ,2 ,1 ,)1(kkktNihixki第8章 有限差分方法 节点处的函数为u(xi,tk)=。在(i,k)点，用中心差商代替,则式(8.86)中的偏微分方程变为 (8.88),22tukiu22xukikiiikikikifuuuhcuuu)2()2(111221

34、12第8章 有限差分方法 引入P=(c/h)2，上式变为 (8.89)由k-1和k时刻的u可直接求k+1时刻的u，不必解联立方程组，故这种差分方程格式是显式的。误差为O(2)+O(h2)，可以证明，当P1时，这种格式是收敛和稳定的。12111 2(1)()kkkkkkiiiiiiuP uP uufu第8章 有限差分方法 8.8.2 初值、边界条件的差分格式初值、边界条件的差分格式初值条件用前向差商代替，有 (8.90)初始条件变为 (8.91)txu)0,()(1)0,(01iiuutxu iiiiixRxRuxRu21110),(第8章 有限差分方法 边界条件可以参照8.4节中的处理方法，类

35、似于式(8.39)可得边界条件的差分式 (8.92)22221221131211132423242bhaububhcubhaububhcunnn第8章 有限差分方法 8.8.3 计算程序流程计算程序流程计算程序流程为 dimension x(N),u(N),V(N),W(N),R1(N),R2(N)输入c,a0,h,tmax,(a,b,c)1,2，P=(c/h)2，N=1+a0/h，kmax=tmax/，22。第8章 有限差分方法 结点坐标与初值do 10 i=1,Nx(i)=(i-1)*h,W(i)=R1(xi)(k=0)10 u(i)=R1(xi)+*R2(xi)(k=1)t=输出t,x,

36、udo 40 k=2,kmaxt=k*，t1=t-do 20 i=2,N-120 V(i)=2*(1-P)*u(i)+P*(u(i+1)+u(i-1)+2*f(x(i),t1)-W(i)第8章 有限差分方法 通过用V表示的边界条件(8.92)式给出V(1)和V(N)。do 30 i=1,NW(i)=u(i)(为下一步准备)30 u(i)=V(i)40 输出t,u对于二维和三维波动方程，可采用交替方向的方法来处理。第8章 有限差分方法【例例8.4】计算一维波动方程2222,01,0(,0)sin0;1(,0)(1)(0,)(1,)0,0yyxttxy xxxl ly xxxtytytt第8章 有

37、限差分方法 对于一维波动方程，可采取如下差分格式由，将上式代入并整理后，得112222kkkiiiyyyyt其中,。222221yyvtx211222kkkiiiyyyyxh222221yyvtx1221112(1)()kkkkkiiiiiyyyyyvh第8章 有限差分方法(1)取=1,h=0.2，计算k=2,3,4,5层的值；(2)取=1,h=0.05，分别用两种差分格式计算k=2,3,20层的近似值。解解(1)由=1,h=0.2及v=1得。由l=1，得。按差分格式有0.2hv5lNh11111216,2,3,4,52,3,4sin(1),2,3,4,5sin(1)(1)(1(1),2,3,

38、4,50kkkkiiiiiikkyyyyikyihiyihihihiyy第8章 有限差分方法 当k=1时，有 即52sin1,401,1yy54sin5sin1,51,2yysin52sin1,61,3yy9511.001411yy5878.05878.01512yy09511.01613yy第8章 有限差分方法 当k=2时，有当k=3时，有9991.002421yy6198.06198.02522yy09991.02623yy6678.003431yy4113.04113.03532yy06678.03633yy第8章 有限差分方法 当k=4时，有当k=5时，有6678.004441yy04

39、80.00480.04542yy06678.04643yy5398.005451yy3313.03313.05552yy05598.05653yy第8章 有限差分方法 计算程序如下:real u(6,5),a,h,v,tparameter pi=3.1415926a=1.h=0.2v=1.tao=a*h/v do i=2,5 第8章 有限差分方法 u(i,1)=sin(i-1)*pi*h)u(i,2)=sin(i-1)*pi*h)+(i-1)*h*tao*(1-(i-1)*h)end dodo j=1,5 u(1,j)=0.u(6,j)=0.end dodo j=2,4do i=2,5 u(i

40、,j+1)=u(i+1,j)+u(i-1,j)-u(i,j-1)第8章 有限差分方法 end doend dodo j=1,5do i=1,6write(*,*)u(i,j)end doend doend 第8章 有限差分方法 程序中a、v、tao分别对应于、v和。第(2)问在此不作解答，请读者自行推导,并参照第(1)问源程序编程。第8章 有限差分方法 8.1 试采用有限差分法编程求解例8.2并统计迭代次数，其中取h=0.1，=1.25，=10-4。习习 题题 八八 第8章 有限差分方法 8.2 试采用有限差分法编程求解泊松方程在网格结点处的值取h=0.2,=1.25，=10-4,并与解析解u

41、=xey作比较。xyyyxxyuxuuuyxxyuxue|,|e2|,0|10,20 e10202222第8章 有限差分方法 8.3 编程完成例8.3中同轴线电势沿径向的分布及特征阻抗计算。8.4 试采用有限差分法编程求解单位圆域中泊松方程在网格结点处的值其中表示半径为1的单位圆边界上u值为1。取h=0.1，=1.25，=10-5，M=16，并与解析解作比较。0),1()2sin(5022uru0),1(u)2sin()1(62522rru第8章 有限差分方法 8.5 编程求解一维扩散方程的解02221110max022 0 e2|,0,0 axcnubuaxcnubuatxuttaxxuDt

42、uxt第8章 有限差分方法 取a1=1,b1=1,c1=0，a2=1,b2=-1，c2=0，a0=1.0，tmax=10，D=0.1，h=0.1,=0.1。输出t=1,2,10时刻的x和u(x)，并与解析解u=exp(x+0.1t)作比较。第8章 有限差分方法 8.6 用ADI格式编程求解二维扩散方程 取h=0.01,=0.05,输出t=1时的u(x,y)结果并与解析解u(x,y,t)=ex+y+t作比较。222200.10.700.1(),00.1,00.7,01(,0)0,0,x yxxx tx tyxuuuDxyttxyu x yeuuuunnuuueuenn 第8章 有限差分方法 8.7 编程计算例8.4题，取=1,h=0.05，计算k=2,3,20的近似值。