微分方程数值解报告

2011-12-22

山东大学数学学院08级基地班 信息与计算科学专业 乔珂欣 学号 200800090114

微分方程数值解报告

目 录

一维变系数二点边值问题的中心差分数值解法 ....................................................... 3

1.中心差分格式的建立 ....................................................................................................................... 3 2.算例 ................................................................................................................................................... 5

二维常系数椭圆问题五点中心差分 ........................................................................... 7

1.建立五点差分格式 ........................................................................................................................... 7 Jacobi迭代 ...................................................................................................................................... 7 SOR迭代 ......................................................................................................................................... 8 GS迭代 ............................................................................................................................................ 8 2.算例及三种迭代格式的性能比较 ................................................................................................... 8 3.交替方向迭代法 ............................................................................................................................. 10 1.交替方向法格式 ......................................................................................................................... 10 2.算例 ............................................................................................................................................. 11 不同迭代法比较 ................................................................................................................................ 13

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

1.

一维抛物方程的初边值问题 .................................................................................................... 14

向后欧拉格式 ................................................................................................................................ 14 算例求解 ........................................................................................................................................ 14 2. 二维抛物方程的初边值问题 .................................................................................................... 17 向后欧拉格式 ................................................................................................................................ 17 交替方向法 .................................................................................................................................... 17 算例求解 ........................................................................................................................................ 18

两点边值问题有限元法 ............................................................................................. 21

1.一般两点边值问题有限元方程格式 .......................................................................................... 21

从Ritz有限元方法出发 ................................................................................................................ 21 从Galerkin有限元法出发 ............................................................................................................ 22 2.算例求解 ........................................................................................................................................ 22 3.收敛性和误差估计 ....................................................................................................................... 24

二阶椭圆问题有限元法 ............................................................................................. 25

1.POISSON方程三角形单元格式 ...................................................................................................... 25 2.算例求解 ......................................................................................................................................... 27

微分方程数值解报告 乔珂欣 200800090114 2011-12-22

附件：MATLAB程序代码： ................................................................................... 29

代码1 ：一维变系数二点边值问题中心差分格式 ...................................................................... 29 代码2 ：二维常系数椭圆问题五点中心差分格式及迭代法解差分方程组（JACOBI) .............. 30 代码3 ：二维常系数椭圆问题五点中心差分格式及迭代法解差分方程组（SOR） ............... 32 代码4：二维常系数椭圆问题五点中心差分格式及迭代法解差分方程组（G_S） ................... 32 代码5：交替方向迭代法 ................................................................................................................. 32 代码6：一维抛物方程向后差分法 ................................................................................................. 35 代码7：二维抛物方程向后差分法 ................................................................................................. 36 代码8：二维抛物方程交替方向法 ................................................................................................. 38 代码9：RITZ有限元法求解一维问题： ......................................................................................... 40 代码10：GALERKIN有限元法求解一维问题： ............................................................................. 41 代码11：有限元法解二维椭圆方程（三角剖分） ....................................................................... 42

2

微分方程数值解报告 乔珂欣 200800090114 2011-12-22

一维变系数二点边值问题的中心差分数值解法

考虑两点边值问题： (1.1) Lu??ddudu(p(x))?r(x)?q(x)u?f(x),a?x?b dxdxdx(1.2) u(a)??,u(b)??

假定p?C1[a,b],p(x)?pmin?0,r,q,f?C[a,b],?,?是给定的常数。

1.中心差分格式的建立

首先取N?1个节点：

a?x0?x1???xi???xN?b,

将区间I?[a,b]分成N个小区间：

Ii:xi?1?x?xi,i?1,2,?N.

于是得到区间I的一个网格剖分。记hi?xi?xi?1。用Ih表示网格内点x1，x2，?，xN?1的集合，Ih表示内点和界点x0?a,xN?b的集合。 取相邻节点xi?1,xi的中点x?i?121?(xi?1?xi)(i?1,2,?,N)，称为半整数点。由节点 2N?12a?x0?x1?x3???x22i?12??x?xN?b

又构成[a,b]的一个对偶剖分。

用差商代替微商，将方程（1.1）在内点xi离散化. 逼近边值问题（1.1）（1.2）的差分方程为：

Lhu(xi)?? (1.3)

u(xi?1)?u(xi)u(xi)?u(xi?1)2[p1?p1]?i?hi?hi?1i?2hi?1hi2ri[u(xi?1)?u(xi?1)]?qiu(xi)?fihi?hi?1i=1,2,…,N-1，

3

微分方程数值解报告 乔珂欣 200800090114 2011-12-22

u0??,uN??

当网格均匀，即hi?h(i?1,2,?,N)时差分方程（1.3）简化为

(1.4) Lhui??ui?1?ui?11[pu?(p?p)u?pu]?r?qiui?fi 1i?111i1i?1i2i?i?i?i?h2h2222

这相当于用一阶中心差商，二阶中心差商依次代替（1.1）的一阶微商和二阶微商的结果。这个方程就是中心差分格式。 式（1.4）用方程组展开：

(1.5)

11111??（p+rh)u?[(p?p)?q]u?(p?rh110311131)u2?f1222?h22h22h22?????11111??2(p1?rkh)uk?1?[2(p1?p1)?qk]uk?2(p1?rkh)uk?1?fk?hk?22hk?2k?2hk?22 ???????1(p?1rh)u?[1(p?p)?q]u?1(p?1rh)u?f3N?1N?213N?1N?11N?1NN?1222?N?N?N?N?h2hh2?2222这是一个以u1,u2,……，uN?1为未知量的线性方程组。方程组两边同时乘以h，再写成矩阵的形式为：

(1.6)

24

微分方程数值解报告 乔珂欣 200800090114 2011-12-22

??12p?p?qh?p?rh0??001?311?322?22????u?112?p?rhp?p?qh?p?rh??00???1?2255332???u2?222222????????????u11??N?2?0??0?p5?rN?2hp3?p5?qN?2h2?p3?rN?2h????u?N?2N?N?N?22222???N?1???120??00?p?rhp?p?qh?N?1N?1???N?32N?1N?3222??1?2?hf?(p?rh)?111??22??2??hf2??????2??hfN?2??1?h2f?(p?rh)???N?1N?12N?1??2? 解此N-1维的线性方程组，可得数值解。

2.算例

求解如下二点边值问题：

?2d2udu1?2x?u??3x2?4x?3?x2(1.7) ?dxdxx?u(1)?6,u(3)?54?该边值问题的精确解为u?x?2x?3x

32x?[1,3]

取步长h=1/10,1/50,1/100，差分解与精确解的最大绝对误差的数量级依次如下表所示

1/10 1/100 1/250 1/500 步长 最大绝对误差数量级 10?4 10?5 10?6 10?7 由此可见，随着步长的减小，差分解的精度越来越高。

给出当N=10时的数值解与精确解的对比图和误差分布图：

5

微分方程数值解报告 乔珂欣 200800090114 2011-12-22

6

微分方程数值解报告 乔珂欣 200800090114 2011-12-22

二维常系数椭圆问题五点中心差分

考虑二维模型问题：

??u?f(x,y),?? (2.1) ?(x,y)?G?(0,1)?(0,1),?u?0,(x,y)??G?

1.建立五点差分格式

取定沿x轴和y轴的步长h1和h2。做两族与坐标轴平行的直线：

x?ih1,y?jh2,i?0,?1,?,j?0,?1,?.

两族直线的交点记为(xi,yj)或(i,j)。

沿x,y方向分别用二阶中心差商代替uxx,uyy，则得 (2.2) ??huij????ui?1,j?2ui,j?ui?1,jui,j?1?2ui,j?ui,j?1????fij 22h1h2??式中uij表示节点(i,j)上的网函数值。

特别取正方形网格：h1?h2?h，则差分方程（2.2）简化为

1h2u?(u?ui,j?1?ui?1,j?ui,j?1)?fij(2.3) ij4i?1,j 4ui0?u0j?uiN?uNj?0,i,j?1,2,?,N?1Jacobi迭代

在Richardson同步迭代

u(k?1)?u(k)??(Au(k)?b),中取?opt?k?0,1,?,

1，得解（2.3）相应的Jacobi迭代 47

微分方程数值解报告 乔珂欣 200800090114 2011-12-22

1(k?1)(k)k)(k)(k)(k)(k)2uij?uij?(?ui(??u?4u?u?u?hfij)1,ji,j?1iji?1,ji,j?14(2.4)

1k)(k)(k)(k)2?(ui(??u?u?u?hfij),k?0,1,?1,ji,j?1i?1,ji,j?14SOR迭代

SOR迭代格式为

(2.5) uij(k?1)(k)?uij??4(k)k?1)(k?1)(k)(k)2(4uij?ui(?1,j?ui,j?1?ui?1,j?ui,j?1?hfij),k?0,1,?

其中?的取值范围为0???2。最佳松弛因子?opt满足?opt?1。

GS迭代

在迭代式（2.5）中取??1，则得Gauss-Seidel迭代格式。

2.算例及三种迭代格式的性能比较

求解边值问题

???u??2?2e?(x?y)(sin?xcos?y?cos?xsin?y),?(x,y)?G?(0,1)?(0,1),(2.6) ??u?0,(x,y)??G?其精确解为u(x,y)?e?(x?y)

sin?xsin?y

取步长h=1/10,1/20,1/30,1/40,1/50当两次迭代值重合到小数点后六位时停止迭代，用

Jacobi迭代，Gauss-Seidel迭代和SOR迭代（松弛因子取为1.75）求解差分方程，三种迭代法迭代次数与精度的比较如下表所示 步长 1/10 1/20 1/30 1/40 Jacobi 迭代次数 120 371 678 1036 最大误差 1.164 0.294 0.138 0.098 SOR 迭代次数 48 48 110 187 最大误差 1.164 0.293 0.132 0.078 G_S 迭代次数 125 393 741 1134 最大误差 1.164 0.294 0.138 0.098 由上表可见，随着区域剖分的加细，数值解的精度越来越大。下面以步长为1/50，采用SOR迭代（松弛因子取1.75）为例给出数值解、精确解和误差绝对值的网格图：

8

微分方程数值解报告 乔珂欣 200800090114 2011-12-22

9

微分方程数值解报告 乔珂欣 200800090114 2011-12-22

3.交替方向迭代法 1.交替方向法格式

考虑边值问题

-?u?f(x,y),u?0,when的五点差分格式

(2.7) ??huij??(0?x,y?1,x?0,1ory?0,1

ui?1,j?2uij?ui?1,jh2?ui,j?1?2uij?ui,j?1h2)?fij

先定义矩阵L1和L2。对向量u?uij，令

??(L1u)ij?(?ui?1,j?2uij?ui?1,j)/h2,(L2u)ij?(?ui,j?1?2uij?ui,j?1)/h2, (Lhu)ij?(L1u)ij?(L2u)ij,则可将（2.7）写成

Lhu?L1u?L2u?f

10

微分方程数值解报告 乔珂欣 200800090114 2011-12-22

其中f?fij。于是PR迭代为：

??(2.7.3)1(2.7.3)2u1(k?)2?u(k)??k(L1u1(k?)2?L2u(k)?f),?L2u(k?1)?f),k?0,1,?.u(k?1)?u1(k?)2??k(L1u1(k?)2

其中?k是迭代参数。按层合并，得

?(2.7.3)1(2.7.3)?2由u(k)(I??kL1)u1(k?)2?(I??kL2)u(k)??kf,1(k?)2

(I??kL2)u(k?1)?(I??kL1)u??kf,到u(k?1)只需交替解两个具三对角系数阵的方程(2.7.3)1?和(2.7.3)2?，可用消去法解。

2.算例

求解边值问题

???u??2?2e?(x?y)(sin?xcos?y?cos?xsin?y),?(x,y)?G?(0,1)?(0,1),??u?0,(x,y)??G?其精确解为u(x,y)?e?(x?y)sin?xsin?y 在迭代中，最优迭代参数取为?opt?

12h(sin?h)?1，当两次迭代值间重合到小数点后2六位时迭代停止。

下面是取不同步长时PR迭代的迭代次数和与精确解相比最大误差的对比表： 步长 最大误差 迭代次数 1/16 0.459 42 1/32 0.114 80 1/64 0.029 153 1/128 0.00717 290 由上表可见，迭代误差随步长的减小而逐渐减小。 下面给出当步长取为1/60时精确解、数值解和二者误差的图形：

11

微分方程数值解报告 乔珂欣 200800090114 2011-12-22

12

微分方程数值解报告 乔珂欣 200800090114 2011-12-22

不同迭代法比较

下面是交替方向迭代法和Jacobi迭代、SOR迭代的迭代次数和与精确解相比最大误差的对比表： 步长 1/10 1/20 1/30 1/40 由上表可以看出，交替方向迭代法具有更高的收敛速度，且随着步长的减小数值解与精确解的误差逐渐减小。

Jacobi 迭代次数 最大误差 SOR 迭代次数 最大误差 PR 迭代次数 最大误差 120 371 678 1036 1.164 0.294 0.138 0.098 48 48 110 187 1.164 0.293 0.132 0.078 27 52 76 99 1.164 0.293 0.13 0.073 13

微分方程数值解报告 乔珂欣 200800090114 2011-12-22

抛物型方程的有限差分法

1. 一维抛物方程的初边值问题

一维热传导方程：

(3.1)

?u?2u?a2?f(x),?t?x0?t?T,

其中a是常数，f(x)是给定的连续函数。按照初值条件的不同给发，可将(3.1)的定解问题分

为两类：

第一， 初值问题：

(3.2)

第二， 初边值问题：

u(x,0)??(x),???x??

(3.3) (3.4)

u(x,0)??(x),u(0,t)?u(l,t)?0,0?x?l 0?t?T

向后欧拉格式

向后差分格式为

(3.5) (3.6)

?1un?unjj??a?1n?1?1un?unj?1?2ujj?1h2?fj

u0j??j??(xj),nnu0?uJ?0,

其中j?1,2,?,J?1,n?0,1,?,N?1.将(3.5)改写为

(3.7)

此差分格式为隐格式。

?1n?1?1n?run?runj?1?(1?2r)ujj?1?uj??fi

算例求解

求解一维抛物方程的初边值问题：

(3.8)

?u?2u??sint,?t?x214

0?x?1,t?0

微分方程数值解报告 乔珂欣 200800090114 2011-12-22

(3.9) (3.10)

ux(0,t)?ux(1,t)?0,u(x,0)?cos?x,t?0, 0?x?1

2设空间步长h=1/J,时间步长??0，tn?n?，网比r??/h.

向后格式为 (3.11)

?1un?unjj???1n?1?1un?unj?1?2ujj?1h2?sintn?1,j?1,?,J?1,n?1,?.N.

边值条件： (3.12)

j?0,n?1nu0?u0?n?1nuJ?uJn?1n?1u1n?1?2u0?u?1??sintn?1,2hn?1n?1n?1uJ?uJ?1?2uJ?1??sintn?1,2hn?1n?1u?1?u1,

(3.13) 初值条件： (3.14)

j?J,?n?1n?1uJ?1?uJ?1,

u(xj,0)?cos?xj,j?0,1,?,J.

计算方案：

A） h=1/40,?=1/1600,此时r=1，计算到时间层t1600?1， B） h=1/80,??1/3200,此时r=1，计算到时间层t3200?1。

两种方案的误差阶都为1?10。以方案A为例，给出精确解与差分解的对比图及误差分布图如下：

?315

微分方程数值解报告 乔珂欣 200800090114 2011-12-22

Matlab程序见附件

16

微分方程数值解报告 乔珂欣 200800090114 2011-12-22

2. 二维抛物方程的初边值问题

考虑平面有界域G上的抛物方程：

(3.15)

?u??(a?u)?f,?t(x,y)?G,0?t?T,

其中a,f都是区域G上给定的光滑函数，a?a(x,y)???0.初值条件为

(3.16)

边值条件为下列三种类型之一： (3.17) (3.18)

u(x,y,0)??(x,y),

u?G??(x,y), （第一边值条件） ?u????(x,y), （第二边值条件）

?G(3.19)

?u??u?G??(x,y), （第三边值条件） ??其中?是边界?G的单位外法向。

向后欧拉格式

(3.20)

2r1??/h12,r2??/h2.

交替方向法

考虑二维热传导方程的边值问题：

(3.21)

?n?11n?1n?1n?1?uij?r1?a1(uin???u)?a(u?u1,ji,j1i,ji?1,j)??i?,ji?,j?22??1n?1n?1n?1?nr2?a1(uin,??u)?a(u?uj?1i,j1i,ji,j?1)??ui,j??fij,i,j??i,j?22?ut?uxx?uyy,0?x,y?l,t?0,??u(x,y,0)??(x,y), ??u(0,y,t)?u(l,y,t)?u(x,0,t)?u(x,l,t)?0.?

其中

取空间步长h=1/M,时间步长??0，作两族平行于坐标轴的网线。由地n层到第n+1层计算分成两步：先由第n层到第n?1层，对uxx用向后差分逼近，对uyy用向前差分逼近，217

微分方程数值解报告 乔珂欣 200800090114 2011-12-22

然后由第n?1层到n+1层，对uxx用向前差分逼近，对uyy用向后差分逼近，得到如下ADI2格式（PR格式）：

(3.22) (3.23)

ujk?unjkn?12?/2un?1jk12n?1?2(?xujk2??y2unjk), h12n?1?1?2(?xujk2??y2unjk), h?ujkn?12?/2ujk?unjkn?12更简单的LOD格式（局部一维格式）：

(3.24) (3.25)

?un?1jk12n?1?2?x(ujk2?unjk), 2h12n?1n?1?2?y(ujk?ujk2), 2h?ujkn?12?算例求解

考虑二维抛物方程初边值问题：

(3.26) (3.27) (3.28) (3.29)

?u?4?2(uxx?uyy),?t(x,y)?G?(0,1)?(0,1),t?0,

u(0,y,t)?u(1,y,t)?0,0?y?1,t?0,

uy(x,0,t)?uy(x,1,t)?0,0?x?1,t?0

u(x,y,0)?sin?xcos?y.

2取空间步长h?h1?h2?1/40，时间步长??1/10，网比r??/h?1，采用向后欧拉格

式和交替方向法（LOD法）分别计算到时间层t=1.下表列出在节点

的计算结果。 方法 向后欧拉法 LOD法 精确解 向后欧拉法 LOD法 精确解 jk(xj,yk)?(,),j,k?1,2,3.

44(xj,yk) y1 x1 0.152931 0.189727 0.145606 0.000000 0.000000 0.000000 x2 0.216277 0.268315 0.205919 0.000000 0.000000 0.000000 18

x3 0.152931 0.189727 0.145606 0.000000 0.000000 0.000000 误差阶 0.01 0.01 0.01 0.01 y2

微分方程数值解报告 乔珂欣 200800090114 2011-12-22

向后欧拉法 LOD法 精确解 y3 -0.152931 -0.189727 -0.145606 -0.216278 -0.268315 -0.205919 -0.152931 -0.189727 -0.145606 0.01 0.01 下面给出该问题精确解、用向后欧拉法所得数值解及用LOD法所得数值解的效果图

19

微分方程数值解报告 乔珂欣 200800090114 2011-12-22

Matlab程序见附件

20

微分方程数值解报告 乔珂欣 200800090114 2011-12-22

两点边值问题有限元法

1.一般两点边值问题有限元方程格式

考虑一般两点边值问题

(4.1)

其中

ddu??Lu??(p)?qu?f,a?x?b, dxdx??u(a)?0,u?(b)?0,?

p?C1(I),p?p0?0,q?C(I),q?0,f?C(I),I??a,b?。

从Ritz有限元方法出发

单元刚度矩阵为

(4.2) K(i)i)?ai(?1,i?1???a(i)?i,i?1i)?ai(?1,i (i)??aii?(4.3)

?a(i)?1[h?1p(x?h?)?hq(x?h?)(1??)2]d?,i?1,i?1ii?1i?0ii?1i?1?(i)aii??[hi?1p(xi?1?hi?)?hiq(xi?1?hi?)?2]d?, ?0?1i)(i)?1?ai(??a?[?hp(xi?1?hi?)?hiq(xi?1?hi?)?(1??)]d?.1,ii,i?1i?0?按规则组装成总刚度矩阵K? 令

?Ki?1n(i)。

(4.4)

以及b?(b1,b2,?,bn)T,

(i)(i)Tf(i)?(fi?1,fi),

?f(i)?h1f(x?h?)(1??)d?,i?0i?1i?i?1 ?1?fi(i)?hi?f(xi?1?hi?)?d?,0?(4.5)

则有限元方程为

b1?f1(1)?f1(2),b2?f2(2)?f2(3),?,bn?fn(n),

Ku?b.

21

(4.6)

微分方程数值解报告 乔珂欣 200800090114 2011-12-22

从Galerkin有限元法出发

Galerkin有限元方程为

(4.7)

?a(?,?)u??ijii?1nbaf?jdx,j?1,2,?,n.

系数矩阵第j行只有三个非零元素，即

(4.8) (4.9) (4.10)

1a(?j?1,?j)??[?h?jp(xj?1?hj?)?hjq(xj?1?hj?)(1??)?]d?,

012a(?j,?j)??[h?jp(xj?1?hj?)?hjq(xj?1?hj?)?]d??011?[h01?1j?1p(xj?hj?1?)?hj?1q(xj?hj?1?)(1??)]d?,102

1a(?j?1,?j)??[?h?j?1p(xj?hj?1?)?hj?1q(xj?hj?1?)(1??)?]d?,

这里j?2,3,?,n?1.第一行只有两个非零元素：a(?1,?1),a(?1,?2).第n行只有两个非零元素：a(?n,?n?1)和

(4.11) (4.12)

?1a(?n,?n)??[hnp(xn?1?hn?)?hnq(xn?1?hn?)?2]d?

01方程(4.7)的右端项

(f,?j)?hj?f(xj?1?hj?)?d??hj?1?f(xj?hj?1?)(1??)d?0011

2.算例求解

考虑常系数两点边值问题

(4.13)

其精确解为y?sin?y????24y??22sin?2x,0?x?1,y(0)?0,y?(1)?0

?2x。算例中p(x)?1,q(x)??24,f(x)??22sin?2x。

将积分区间等分为N份，则步长hi?刚度矩阵的元素为

1,i?1,2,?n，记为h。从Ritz法出发，单元N(4.14)

1i)?12?ai(??[hp?hq(1??)]d?,1,i?1?0?1?(i)?12a?[hp?hq?]d?, ?ii?0?1i)(i)?1?ai(?1,i?ai,i?1??0[?hp?hq?(1??)]d?.?22

微分方程数值解报告 乔珂欣 200800090114 2011-12-22

对应积分式(4.4)的积分式为

21??(i)?f?hsin[(xi?1?h?)](1??)d?,i?1??0?22 ?2?f(i)?h1?sin[?(x?h?)]?d?,i?022i?1??(4.15)

积分采用matlab自带的积分函数quad。步长依次取为1/5,1/10,1/50,1/100， 最大误差按

精确解与数值解之差的绝对值的最大值记，得如下结果：

表格 1 步长(h) 最大误差(err)

1/5 0.004109 1/10 0.001028 1/50 0.000041 1/100 0.000010 以步长取为h=1/50为例，给出从Ritz法出发的有限元法数值解与精确解的对比图如下：

图 0-1

从Galerkin法出发，对应(4.8)(4.9)(4.10)的系数矩阵元素依次为

(4.16) (4.17) (4.18)

a(?j?1,?j)??[?h?1p?hq(1??)?]d?,

01a(?j,?j)??[hp?hq?]d???[h?1p?hq(1??)2]d?,

?120011a(?j?1,?j)??[?h?1p?hq(1??)?]d?,

01j?2,3,?,n?1.对应(4.11)和(4.12)的元素分别为

23

微分方程数值解报告 乔珂欣 200800090114 2011-12-22

(4.19)

(4.20)(f,?j)?h?01a(?n,?n)??[h?1p?hq?2]d?01

??2sin[(xj?1?h?)]?d??h22??1?20sin[(xj?h?)](1??)d?

22积分采用matlab自带的积分函数quad.步长依次取为1/5,1/10,1/50,1/100,1/200， 最大误

差按精确解与数值解之差的绝对值的最大值记，得如下结果：

表格 2

步长（h） 1/5 1/10 0.144690 1/50 0.028851 1/100 0.014416 1/200 0.07206 最大误差（err） 0.289031 以h=1/50为例给出Galerkin有限元法数值解和精确解的图形如下图：

图 0-2

（相关程序见附录）

3.收敛性和误差估计

令u和uh分别表示精确解和有限元解在剖分区间节点处的值，收敛性表示为

u(xi)?uh(xi) (4.21) Mhk?u?uh?0max?i?N?1记最大误差为err,则问题转化为在方程

(4.22)

Mhk?err

24

微分方程数值解报告 乔珂欣 200800090114 2011-12-22

中，已知h和err，求解M和k的拟合问题。在Matlab中拟合采用最小二乘法实现。

取表1中的四组数据，对err和h进行最小二乘幂函数拟合，求得从Ritz法出发的误差 阶为k=2.019,M=0.106.

取表2中的五组数据，对err和h进行最小二乘幂函数拟合，求得从Galerkin法出发的误差阶为k=1.001,M=1.447.

二阶椭圆问题有限元法

1. Poisson方程三角形单元格式

考虑Poisson 方程第一边值问题

(5.1) ??u?f(x,y),?????(x,y)?G,

G是xy平面的一个有界域，其边界?是分段光滑的简单闭曲线。在?上给出第一边值条件

(5.2)

u????(x,y).

其中，f(x,y),?(x,y)是给定的连续函数.

以Poisson 方程第一边值问题(5.1),(5.2)为例说明用有限元法解题的主要过程. 第一， 把边值问题(5.1),(5.2)化成等价的变分形式：求u?H1(G),u??=?使

(5.3)

其中， (5.4)

1a(u,v)?f(u,v),对任意v?H0(G),

??u?v?u?v?a(u,v)?????dxdy. ??x?x?y?y?G?第二，对区域G作网格剖分.这里我们用三角剖分. 首先在?上取有限个点，依次联成一闭多边形?h(h表示单元的最大直径)，以之近似代替?, 并以?h围成的多边形域Gh近似代替G.然后，把Gh分割成有限个三角单元之和, 从而得到Gh的一个三角剖分. 确定 剖分后，按一定次序单元顶点编号. 第i号点的坐标记为(xi,yi). 设?=(?1,?2,?3)，用

??=?(?1,?2,?3)表示一个三角元，其顶点编号为?1,?2,?3，并且顶点排列顺序?1,?,?3

符合逆时针方向. S?表示??的面积.

第三，构造奇函数或单元形状函数. 我们采取线性元，则每一节点i对应一基函数?i.

?i的

25

微分方程数值解报告 乔珂欣 200800090114 2011-12-22

支集是一切以i为顶点的三角单元 ??=?(i,j,k)(?1?i,?2?j,?3?k). 为

?i在??的限制

(5.5)

其中，

?i?Li?1(ai?bix?ciy), 2S?(5.6) (5.7)

ai?xjyk?xkyj,???bi?yj?yk,???ci?xk?xj,

1xi2S??1xj1xkyiyj. yk假若与i相邻的所有节点是le(e?1,2,?，m), 则只要（5.5）~（5.7）中依次取j?le,?k?le?1， 其中e?1,2,?m,??lm?1?l1 （i是内点），或e?1,2,?,m?1（i是界点），就得出?i在支集上的表达式.

第四, 形成有限元方程. 假定内点的编号是1,2，?,n1, 界点的编号是

n1?1,n2?1,?,n1?n2, 则有限元方程为

(5.8) 其中 (5.9)

?a(??)ui,ji?1n1i?(f,?j)?n1?n2i?n1?1?a(?,?)?,???j?1,2?,n,

iji1?i??(xi?yi),??i?n1?1,?,n1?n2,

(f,?j)???f?jdxdy?G11?2?S???(a??j?bjx?cjy)fdxdy.

??表示就所有以j为顶点的??求和, 而

(5.10)

???i??j??i??j?a(?i,?j)????????dxdy?x?x?y?y?G?1 ?????????????????2??(bibj?cicj)dxdy4S????????????????????11?(bibj?cicj),4?S?在程序编写时，先形成单元刚度矩阵，然后组合成总刚度矩阵。Matlab程序见附件。

26

微分方程数值解报告 乔珂欣 200800090114 2011-12-22

2.算例求解

???u??2?2e?(x?y)(sin?xcos?y?cos?xsin?y),?(x,y)?G?(0,1)?(0,1),??u?0,(x,y)??G?(5.11)

其精确解为u(x,y)?e?(x?y)sin?xsin?y。

当剖分为312个单元时，数值解与精确解的最大绝对值误差为0.840078；当剖分加细，

剖分为1248个单元时，误差为0.328238.下面为剖分单元为312时数值解、精确解和误差分布图。

27

微分方程数值解报告 乔珂欣 200800090114 2011-12-22

28

微分方程数值解报告 乔珂欣 200800090114 2011-12-22

附件：Matlab程序代码：

代码1 ：一维变系数二点边值问题中心差分格式

%一维变系数二点边值问题中心差分格式 clc clf syms x;

a=1; %区间界点 %a=0;

b=3; %区间界点 %b=1

p=-x^2; %这是p函数 %p=exp(x); q=1/x; %这是q函数 %q=sin(x)+1+x;

f=-3*x^2-4*x+3;%这是f函数

%f=-exp(x)*(2*x+1)+(sin(x)+1+x)*x*(x-1); r=-4*x; %这是r函数. %r=0;

N=10 %将区间划分的等分,这里控制！ h=(b-a)/N; %这里确定步长 value_of_f=zeros(N-1,1);%这是f diag_0=zeros(N-1,1);%确定A的对角元

diag_1=zeros(N-2,1);%确定A的偏离对角的上对角元 diag_2=zeros(N-2,1);%确定A的偏离对角的下对角元 X=a:h:b;

u_a=6; %边界条件 %u_a=0;

u_b=54; %边界条件 %u_b=2; for j=1:N-1

diag_0(j)=(subs(p,{x},{(X(j)+X(j+1))/2}))+(subs(p,{x},{(X(j)+X(j+1))/2+h}))+(subs(q,{x},{X(j+1)}))*(h^2);

end %获取对角元素 for j=1:N-2

diag_2(j)=-(subs(p,{x},{(X(j+1)+X(j+2))/2}))-((subs(r,{x},{X(j+2)}))*h)/2; end %获取A的第三条对角 for j=1:N-2

diag_1(j)=-(subs(p,{x},{(X(j+1)+X(j+2))/2}))+((subs(r,{x},{X(j+1)}))*h)/2; end %获取A的第二条对角 for j=1:N-1;

29

微分方程数值解报告 乔珂欣 200800090114 2011-12-22

value_of_f(j)=(h^2)*(subs(f,{x},{X(j+1)})); end %获取F值

value_of_f(1)=value_of_f(1)+u_a*(subs(p,{x},{(X(2)+X(1))/2})+(subs(r,{x},{X(2)}))*h/2);

value_of_f(N-1)=value_of_f(N-1)+u_b*(subs(p,{x},{(X(N)+X(N+1))/2})-(subs(r,{x},{X(N)}))*h/2);

A=diag(diag_0)+diag(diag_1,1)+diag(diag_2,-1);%组装系数矩阵 format long e

U=inv(A)*value_of_f ; %差分解 fprintf('\\n'); dx=X(2:N);

precise_value=dx.^3+2.*(dx.^2)+3.*dx; %精确解 %precise_value=dx.*(dx-1);

deta=U-precise_value' ; %误差 deta_max=max(abs(deta));%最大误差 fprintf('最大的误差是%f\\n',deta_max) figure(1);

plot(X(2:N),U,'b*',X(2:N),precise_value,'r--') ; %差分解与精确解对比表 title('差分解与精确解对比图'); legend('差分解','精确解'); figure(2);

plot(X(2:N),deta) ;%误差图 title('误差图');

代码2 ：二维常系数椭圆问题五点中心差分格式及迭代法解差分方程组（Jacobi)

%二维常系数椭圆问题五点中心差分格式及迭代法解差分方程组（Jacobi) clc clf clf syms x; syms y;

x_a=0; x_b=1;%x的区间端点 y_a=0;y_b=1;%y的区间端点

f=-2*(pi^2)*exp(pi*(x+y))*(sin(pi*x)*cos(pi*y)+cos(pi*x)*sin(pi*y));%这是f函数 N=10%将区间划分的等份，这里控制！

h=1/N;%步长。就本题而言，x,y方向的步长是一样的 U=ones(N+1,N+1)*10;%迭代矩阵，含边界点的值 %按题意将边界点的值取为0 for j=1:N+1 U(1,j)=0;

30

微分方程数值解报告 乔珂欣 200800090114 2011-12-22

U(N+1,j)=0; U(j,1)=0; U(j,N+1)=0; end

%构造f的网格函数 X=x_a:h:x_b; Y=y_a:h:y_b;

values_of_f=zeros(N+1,N+1); for i=1:N+1 for j=1:N+1

values_of_f(i,j)=subs(f,{x,y},{X(i),Y(j)}); end end

%Jacobi迭代 w=1/4;%最佳松弛因子 k=1;%用于统计迭代次数 for i=2:N for j=2:N U1=U;

U(i,j)=w*(U1(i-1,j)+U1(i,j-1)+U1(i+1,j)+U1(i,j+1)+(h^2)*values_of_f(i,j

));

end end

while (max(max(abs(U1-U)))>1e-6)&&(k<1e5) for i=2:N for j=2:N

U1=U;

U(i,j)=w*(U1(i-1,j)+U1(i,j-1)+U1(i+1,j)+U1(i,j+1)+(h^2)*values_of_f(i,j));

end end k=k+1; end

precise_value=exp(pi*(x+y))*sin(pi*x)*sin(pi*y);%精确解 %构造精确解网格函数

precisevalues=zeros(N+1,N+1); for i=2:N for j=2:N

precisevalues(i,j)=subs(precise_value,{x,y},{X(i),Y(j)}); end end

deta=abs(U-precisevalues);%绝对误差

31

微分方程数值解报告 乔珂欣 200800090114 2011-12-22

deta_max=max(max(deta));

fprintf('Jacobi迭代的最大的误差是%f\\n',deta_max) fprintf('Jacobi迭代的迭代次数是%f\\n',k) figure(1);

[x_l,y_l]=meshgrid(X);%绝对误差用 mesh(x_l,y_l,deta); title('误差网格分布'); figure(2);

mesh(x_l,y_l,precisevalues);%精确值的网格函数值 title('精确解'); figure(3);

mesh(x_l,y_l,U);%数值解的网格函数 title('数值解'); U;

代码3 ：二维常系数椭圆问题五点中心差分格式及迭代法解差分方程组（SOR）

上例中w取1.75

代码4：二维常系数椭圆问题五点中心差分格式及迭代法解差分方程组（G_S）

上例中w取1

代码5：交替方向迭代法

%交替方向迭代法 clf clf syms x; syms y;

x_a=0; x_b=1;%x的区间端点 y_a=0;y_b=1;%y的区间端点

f=-2*pi^2*exp(pi*(x+y))*(sin(pi*x)*cos(pi*y)+cos(pi*x)*sin(pi*y));%这是f函数 N=60%将区间划分的等份，这里控制！

h=1/N;%步长。就本题而言，x,y方向的步长是一样的 U=ones(N+1,N+1);%迭代矩阵，含边界点的值 %按题意将边界点的值取为0

32

微分方程数值解报告 乔珂欣 200800090114 2011-12-22

for j=1:N+1 U(1,j)=0; U(N+1,j)=0; U(j,1)=0; U(j,N+1)=0; end

%构造f的网格函数 X=x_a:h:x_b; Y=y_a:h:y_b;

values_of_f=zeros(N+1,N+1); for i=1:N+1 for j=1:N+1

values_of_f(i,j)=subs(f,{x,y},{X(i),Y(j)}); end end

%交替方向迭代

t=h^2/sin(pi*h)/2; %迭代参数 diag_0=zeros(N-1,1); %A的主对角线 diag_1=zeros(N-2,1);%A的次对角线 for i=1:N-1

diag_0(i)=1+2*t/(h^2); end

for i=1:N-2

diag_1(i)=-t/(h^2); end

A=diag(diag_0)+diag(diag_1,1)+diag(diag_1,-1);%组装系数矩阵 U0=U(2:N,2:N); %第一次迭代

%构造第一个迭代方程组的右端向量 f1=values_of_f(2:N,2:N); for i=1:N-1 for j=1:N-1

f1(i,j)=(1-2*t/(h^2))*U(i+1,j+1)+(t/(h^2))*U(i+1,j)+(t/(h^2))*U(i+1,j+2)+t*values_of_f(i+1,j+1);

end end

for j=1:N-1

f1(1,j)=f1(1,j)+t/(h^2)*U(1,j+1); f1(N-1,j)=f1(N-1,j)+t/(h^2)*U(N+1,j+1); end

U1=zeros(N-1,N-1);

for i=1:N-1 %解关于x的方程组 U1(:,i)=inv(A)*f1(:,i);

33

微分方程数值解报告 乔珂欣 200800090114 2011-12-22

end

U(2:N,2:N)=U1;

%构造第二个迭代方程组的右端向量 f2=values_of_f(2:N,2:N); for i=1:N-1 for j=1:N-1

f2(i,j)=(1-2*t/(h^2))*U(i+1,j+1)+(t/(h^2))*U(i,j+1)+(t/(h^2))*U(i+2,j+1)+t*values_of_f(i+1,j+1);

end end

for i=1:N-1

f2(i,1)=f2(i,1)+t/(h^2)*U(i+1,1); f2(i,N-1)=f2(i,N-1)+t/(h^2)*U(i+1,N+1); end

for i=1:N-1 %解关于y的方程组 U1(i,:)=(inv(A)*(f2(i,:)'))'; end k=1;

while(max(max(abs(U1-U0)))>1e-6)&&(k<1e5) U0=U1;

U(2:N,2:N)=U1; for i=1:N-1 for j=1:N-1

f1(i,j)=(1-2*t/(h^2))*U(i+1,j+1)+(t/(h^2))*U(i+1,j)+(t/(h^2))*U(i+1,j+2)+t*values_of_f(i+1,j+1);

end end for j=1:N-1

f1(1,j)=f1(1,j)+t/(h^2)*U(1,j+1); f1(N-1,j)=f1(N-1,j)+t/(h^2)*U(N+1,j+1); end

for i=1:N-1

U1(:,i)=inv(A)*f1(:,i); end

U(2:N,2:N)=U1; %此时得到的是U（k+1/2) f2=values_of_f(2:N,2:N); for i=1:N-1 for j=1:N-1

f2(i,j)=(1-2*t/(h^2))*U(i+1,j+1)+(t/(h^2))*U(i,j+1)+(t/(h^2))*U(i+2,j+1)+t*values_of_f(i+1,j+1);

end

34

微分方程数值解报告 乔珂欣 200800090114 2011-12-22

end

for i=1:N-1

f2(i,1)=f2(i,1)+t/(h^2)*U(i+1,1); f2(i,N-1)=f2(i,N-1)+t/(h^2)*U(i+1,N+1); end

for i=1:N-1

U1(i,:)=(inv(A)*(f2(i,:)'))'; end k=k+1; end

U(2:N,2:N)=U1;

precise_value=exp(pi*(x+y))*sin(pi*x)*sin(pi*y);%精确解 %构造精确解网格函数

precisevalues=zeros(N+1,N+1); for i=2:N for j=2:N

precisevalues(i,j)=subs(precise_value,{x,y},{X(i),Y(j)}); end end

deta=abs(U-precisevalues);%绝对误差 deta_max=max(max(deta));

fprintf('交替方向迭代的最大的误差是%f\\n',deta_max) fprintf('交替方向迭代的迭代次数是%f\\n',k) figure(1);

[x_l,y_l]=meshgrid(X);%绝对误差用 mesh(x_l,y_l,deta); title('误差网格分布'); figure(2);

mesh(x_l,y_l,precisevalues);%精确值的网格函数值 title('精确解'); figure(3);

mesh(x_l,y_l,U);%数值解的网格函数 title('数值解'); U;

代码6：一维抛物方程向后差分法

%一维抛物方程向后差分 J=80;%区间剖分，这里控制！ h=1/J;

X=0:h:1;%将区间[0,1]进行剖分

n=3200;%到时间t=1的迭代次数，这里控制！ t=1/n;%时间步长 r=t/(h^2); %网比

35

微分方程数值解报告 乔珂欣 200800090114 2011-12-22

diag_0=zeros(J+1,1); diag_1=zeros(J,1); for i=1:J+1

diag_0(i)=1+2*r; end for i=1:J

diag_1(i)=-r; end

A=diag(diag_0)+diag(diag_1,-1)+diag(diag_1,1);%组装迭代系数矩阵 A(1,2)=-2*r; %边界条件使然 A(J+1,J)=-2*r; %边界条件使然 U=zeros(J+1,1); for i=1:J+1 %初值条件 U(i)=cos(pi*X(i)); end

b=zeros(J+1,1); for i=1:n

for j=1:J+1 %构造右端向量 b(j)=U(j)+t*sin(t*i); end

U=inv(A)*b; end

precise_value=exp(-pi*pi)*cos(pi.*X)+1-cos(1); deta=U-precise_value' ; %误差 deta_max=max(abs(deta));%最大误差 fprintf('最大的误差是%f\\n',deta_max) figure(1);

plot(X,U,'b*',X,precise_value,'r--') ; %差分解与精确解对比表 %axis equal;

title('差分解与精确解对比图'); legend('差分解','精确解'); figure(2);

plot(X,deta) ;%误差图 title('误差图');

代码7：二维抛物方程向后差分法

clc clear all J=40; h=1/J; X=0:h:1; Y=0:h:1;

36

微分方程数值解报告 乔珂欣 200800090114 2011-12-22

N=10; t=1/N; a=1/16; r=a*t/(h^2); A=zeros((J+1)^2); U=zeros((J+1)^2,1); %构造系数矩阵 B=zeros(J+1); for i=2:J

B(i,i)=1+4*r; B(i,i-1)=-r; B(i,i+1)=-r; end B(1,1)=1; B(J+1,J+1)=1; C=-r*eye(J+1); C(1,1)=0; C(J+1,J+1)=0; for i=J+2:J+1:J*J A(i:i+J,i:i+J)=B;

A(i:i+J,i+J+1:i+2*J+1)=C; A(i:i+J,i-J-1:i-1)=C; end

A(1:J+1,1:J+1)=eye(J+1);

A(J^2+J+1:(J+1)^2,J^2+J+1:(J+1)^2)=eye(J+1); A(1:J+1,J+2:2*J+2)=-eye(J+1);

A(J^2+J+1:(J+1)^2,J^2:J^2+J)=-eye(J+1); A(1,J+2)=0;

A((J+1)^2,J^2+J)=0; %构造初始右端向量 for i=2:J

k=(i-1)*(J+1); for j=2:J

U(k+j)=sin(pi*X(j))*cos(pi*Y(i)); end end %

for m=1:N U1=inv(A)*U; for i=2:J

k=(i-1)*(J+1); for j=2:J U(k+j)=U1(k+j); end

37

微分方程数值解报告 乔珂欣 200800090114 2011-12-22

end end

D=zeros(J+1); for i=1:J+1 k=(i-1)*(J+1); D(i,:)=U1(k+1:k+1+J); end

%构造精确解函数 for i=1:J+1 for j=1:J+1

precise_value(i,j)=sin(pi.*X(j))*cos(pi.*Y(i))*exp(-pi^2/8); end end figure(1);

[X,Y]=meshgrid(X,Y); surf(X,Y,precise_value); title('精确解曲面'); figure(2); surf(X,Y,D); title('数值解');

deta=max(max(abs(D-precise_value))); fprintf('最大误差为%f\\n',deta)

代码8：二维抛物方程交替方向法

%交替方向差分 clc clf

x_a=0; x_b=1;%x的区间端点 y_a=0;y_b=1;%y的区间端点 N=40;%控制空间区域划分 h=1/N;%空间步长 x=[x_a:h:x_b]; y=[y_a:h:y_b]; N2=1600;

tao=1/N2;%时间步长 r=tao/(h^2);%网比 a=1/16;

U=ones(N+1,N+1);%迭代矩阵 %按题意将边界点的值取为0 for j=1:N+1 U(1,j)=0; U(N+1,j)=0;

38

微分方程数值解报告 乔珂欣 200800090114 2011-12-22

end %初值条件 for i=2:N for j=1:N+1

U(i,j)=sin(pi*x(i))*cos(pi*y(j)); end end

%差分格式方程组的系数矩阵 diag_0=(1+r*a)*ones(N-1,1); diag_1=(-r*a/2)*ones(N-2,1)';

A=diag(diag_0)+diag(diag_1,1)+diag(diag_1,-1);%组装系数矩阵 A2=zeros(N+1); A2(2:N,2:N)=A; A2(1,1)=1; A2(N+1,N+1)=1; A2(1,2)=-1; A2(N+1,N)=-1; A2(2,1)=-r*a/2; A2(N,N+1)=-r*a/2; f=zeros(N-1,1); f2=zeros(N+1,1);

for n=1:N2 %计算到时间层t=1 %x方向的迭代 for k=2:N

for j=1:N-1 %边界值为0，不必特殊处理j=1和N-1的情况 f(j)=r/32*(U(j,k)+U(j+2,k))+(1-r/16)*U(j+1,k); end

U(2:N,k)=A\\f; end

%y方向的迭代 for j=2:N for k=2:N

f2(k)=r*a/2*(U(j,k+1)+U(j,k-1))+(1-r*a)*U(j,k); end

U(j,:)=(A2\\f2)'; end end

%构造t=1时精确解网格函数 precisevalues=zeros(N+1,N+1); for i=1:N+1 for j=1:N+1

precisevalues(i,j)=sin(pi*x(i))*cos(pi*y(j))*exp(-pi^2/8); end end

39

微分方程数值解报告 乔珂欣 200800090114 2011-12-22

deta=abs(U-precisevalues);%绝对误差 deta_max=max(max(deta));

fprintf('最大误差%f\\n',deta_max) figure(1);

[x_l,y_l]=meshgrid(x);%绝对误差用 mesh(x_l,y_l,deta); title('误差网格分布'); figure(2);

mesh(x_l,y_l,precisevalues');%精确值的网格函数值 title('精确解'); figure(3);

mesh(x_l,y_l,U');%数值解的网格函数 title('数值解'); U;

代码9：Ritz有限元法求解一维问题：

clc clf p=1; %p q=(pi^2)/4; %q

N=50 %剖分份数，这里控制 h=1/N; X=0:h:1;

A=zeros(N+1);%构造系数矩阵 for i=2:N+1

fun1=@(ks)(p./h+h.*q.*((1-ks).^2)); fun2=@(ks)(p./h+h.*q.*(ks.^2)); fun3=@(ks)(-p./h+h.*q.*ks.*(1-ks)); A(i-1,i-1)=A(i-1,i-1)+quad(fun1,0,1); A(i-1,i)=quad(fun3,0,1); A(i,i-1)=quad(fun3,0,1); A(i,i)=quad(fun2,0,1); end

A=A(2:N+1,2:N+1);

values_f=zeros(N+1,1);%构造右端向量 for i=2:N+1

f1=@(ks)(pi.^2./2.*sin(pi./2.*(X(i-1)+h.*ks)).*(1-ks)); f2=@(ks)(pi.^2./2.*sin(pi./2.*(X(i-1)+h.*ks)).*ks); values_f(i-1)=values_f(i-1)+h.*quad(f1,0,1); values_f(i)=h.*quad(f2,0,1); end

values_f=values_f(2:N+1);

40

微分方程数值解报告 乔珂欣 200800090114 2011-12-22

U=A\\values_f; dx=X;

precise_value=sin(pi/2.*dx);%精确解 figure(1);

plot(X,[0;U],'b*',X,precise_value,'r--'); title('Ritz有限元法解与精确解对比图'); legend('数值解','精确解');

err=max(abs([0;U]-precise_value')); format long

sprintf('Ritz有限元法最大误差%f\\n',err)

代码10：Galerkin有限元法求解一维问题：

clf p=1; %p q=(pi^2)/4; %q

N=50 %剖分份数，这里控制 h=1/N; X=0:h:1; A=zeros(N); for i=3:N

fun3=@(ks)-p./h+h.*q.*ks.*(1-ks);

fun2=@(ks)p./h+h.*q.*(ks.^2)+p./h+h.*q.*((1-ks).^2); fun1=@(ks)-p./h+h.*q.*ks.*(1-ks); A(i-1,i-2)=quad(fun1,0,1); A(i-1,i-1)=quad(fun2,0,1); A(i-1,i)=quad(fun3,0,1); end

fun2=@(ks)p./h+h.*q.*(ks.^2)+p./h+h.*q.*((1-ks).^2); A(1,1)=quad(fun2,0,1);

fun3=@(ks)-p./h+h.*q.*ks.*(1-ks); A(1,2)=quad(fun3,0,1);

fun3=@(ks)-p./h+h.*q.*ks.*(1-ks); A(N,N-1)=quad(fun3,0,1); fun4=@(ks)p./h+h.*q.*(ks.^2); A(N,N)=quad(fun4,0,1); values_f=zeros(N,1); for i=2:N+1

f1=@(ks)pi.^2./2.*sin(pi./2.*(X(i)+h.*ks)).*(1-ks)+pi.^2./2.*sin(pi./2.*(X(i-1)+h.*ks)).*ks;

values_f(i-1)=h.*quad(f1,0,1); end

U=inv(A)*values_f;

41

微分方程数值解报告 乔珂欣 200800090114 2011-12-22

dx=X;

precise_value=sin(pi/2.*dx); figure(1);

plot(X,[0;U],'b--',X,precise_value,'r--'); title('Galerkin有限元法解与精确解对比图'); legend('数值解','精确解');

err=max(abs([0;U]-precise_value')); format long

sprintf('Galerkin有限元法最大误差%f\\n',err)

代码11：有限元法解二维椭圆方程（三角剖分）

syms x; syms y;

fct=-2*pi^2*exp(pi*(x+y))*(sin(pi*x)*cos(pi*y)+cos(pi*x)*sin(pi*y));%右端函数 üt=2; point=p; edge=e; trang=t;

n1=size(point,2);%节点个数 n2=size(trang,2);%单元个数 n3=size(edge,2);%边界节点个数 A=zeros(n1); f=zeros(n1,1); for trg=1:n2 i=trang(1,trg); j=trang(2,trg); k=trang(3,trg);

xi=point(1,i);yi=point(2,i); xj=point(1,j);yj=point(2,j); xk=point(1,k);yk=point(2,k); x1=xi*2/3+xj/6+xk/6; y1=yi*2/3+yj/6+yk/6; x2=xi/6+xj*2/3+xk/6; y2=yi/6+yj*2/3+yk/6; x3=xi/6+xj/6+xk*2/3;

y3=yi/6+yj/6+yk*2/3; ai=xj*yk-xk*yj; bi=yj-yk; ci=xk-xj; aj=xk*yi-xi*yk; bj=yk-yi; cj=xi-xk; ak=xi*yj-xj*yi;

42

微分方程数值解报告 乔珂欣 200800090114 2011-12-22

bk=yi-yj; ck=xj-xi; Jcb=[1 xi yi; 1 xj yj; 1 xk yk]; s=det(Jcb)/2; A1=zeros(3); f1=zeros(3,1); %形成单元刚度矩阵 A1(1,1)=bi*bi+ci*ci; A1(1,2)=bi*bj+ci*cj; A1(1,3)=bi*bk+ci*ck; A1(2,1)=bj*bi+cj*ci; A1(2,2)=bj*bj+cj*cj; A1(2,3)=bj*bk+cj*ck; A1(3,1)=bk*bi+ck*ci; A1(3,2)=bk*bj+ck*cj; A1(3,3)=bk*bk+ck*ck; A1=A1./(4*s); %单元右端向量 f1(1)

=(ai+bi*x1+ci*y1).*subs(fct,{x,y},{x1,y1})/3+(ai+bi*x2+ci*y2).*subs(fct,{x,y},{x2,y2})/3+(ai+bi*x3+ci*y3).*subs(fct,{x,y},{x3,y3})/3;

f1(2)=(aj+bj*x1+cj*y1)*subs(fct,{x,y},{x1,y1})/3+(aj+bj*x2+cj*y2)*subs(fct,{x,y},{x2,y2})/3+(aj+bj*x3+cj*y3)*subs(fct,{x,y},{x3,y3})/3;

f1(3)=(ak+bk*x1+ck*y1)*subs(fct,{x,y},{x1,y1})/3+(ak+bk*x2+ck*y2)*subs(fct,{x,y},{x2,y2})/3+(ak+bk*x3+ck*y3)*subs(fct,{x,y},{x3,y3})/3; f1=f1./2;

%组装成总刚度矩阵 A(i,i)=A(i,i)+A1(1,1); A(i,j)=A(i,j)+A1(1,2); A(i,k)=A(i,k)+A1(1,3); A(j,i)=A(j,i)+A1(2,1); A(j,j)=A(j,j)+A1(2,2); A(j,k)=A(j,k)+A1(2,3); A(k,i)=A(k,i)+A1(3,1); A(k,j)=A(k,j)+A1(3,2); A(k,k)=A(k,k)+A1(3,3); %总右端向量 f(i)=f(i)+f1(1); f(j)=f(j)+f1(2); f(k)=f(k)+f1(3);

43

微分方程数值解报告 乔珂欣 200800090114 2011-12-22

end

%边界条件的设定 for m=1:n3

for i=1:n1

A(edge(1,m),i)=0; end

A(edge(1,m),edge(1,m))=1; f(edge(1,m))=0; end U=A\\f;

precise_value=exp(pi*(x+y))*sin(pi*x)*sin(pi*y);%精确解 pv=zeros(n1,1); for i=1:n1

pv(i)=subs(precise_value,{x,y},{point(1,i),point(2,i)}); end

deta=abs(U-pv); deta_max=max(deta);

fprintf('最大误差%f\\n',deta_max) figure(1);

[xx,yy,zz]=griddata(point(1,:)',point(2,:)',U,linspace(min(point(1,:)'),max(point(2,:)'))',linspace(min(point(2,:)'),max(point(2,:)')),'v4'); mesh(xx,yy,zz)

title('有限元数值解'); figure(2);

[xx,yy,zz]=griddata(point(1,:)',point(2,:)',deta,linspace(min(point(1,:)'),max(point(2,:)'))',linspace(min(point(2,:)'),max(point(2,:)')),'v4'); mesh(xx,yy,zz) title('误差分布') figure(3);

[xx,yy,zz]=griddata(point(1,:)',point(2,:)',pv,linspace(min(point(1,:)'),max(point(2,:)'))',linspace(min(point(2,:)'),max(point(2,:)')),'v4'); mesh(xx,yy,zz) title('精确解');

44

本文来源：https://www.bwwdw.com/article/vm2p.html