二阶椭圆偏微分方程实例求解(附matlab代码)

  1. 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
  2. 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
  3. 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。
3.3 编程求解过程中的问题 3.3.1 问题产生 当按照老师要求,n=100,h=1/100 时,运行编好的 matlab 程序时,会出现 如图 3.1 的错误提示。
图 3.1 3.3.2 问题分析 在 matlab 的命令窗口输入“memory” ,出现如图 3.2 的内存使用情况,可以 得出: Memory used by MATLAB: 454 MB (4.760e+008 bytes)。 , 若不用稀疏矩阵定 义 P, 经过粗略计算, 我发现矩阵 P 就要占 800MB 左右的内存, 加上其他数据, 内存消耗至少在 1G 以上。但是我电脑上分配给 matlab 的内存只有:454 MB, 即使在关闭杀毒软件等大部分应用程序后,分配给 matlab 的内存也刚够 1G。
3.2matlab 编程 因为矩阵阶数过大, 所以此题的编程难点为构造系数矩阵,即对线性方程组 的赋值。我采用的方法是分块赋值。对于 P 的赋值,过程如下: 第一步: ������������1 −������������2 bcdi = 0
−������������1 ������������2 ⋱ 0
二、问题分析与模型建立
2.1 微分方程上的符号说明 ������ ������, ������ = ������ ������ ������ ������, ������ = ������ ������ ������ ������, ������ = ������ + ������������ ������, ������ = ������ − ������ ������ ������, ������ = 1������ ������, ������ = y 2 e y x 2 e x e xy y 2 x 2 1 e xy 2.2 课本上差分方程的缺陷 课本上的差分方程为: ������������������ ������������������ − ������������−1,������ ������������−1,������ + ������������ ,������ −1 ������������ ,������ −1 + ������������ +1,������ ������������ +1,������ + ������������ ,������ +1 ������������ ,������ +1 = ������������������
而 f 是(n − 1) 维的列向量,具体如下: ������11 ������12 ⋮ ������ = ������1,������−1 ������21 ⋮ ⋮ ������������ −1,������−1
2
三、求解过程
3.1 对系数矩阵的分析 对上述模型的求解就是对线性方程组的求解。通过观察,我发现 P 是一个 对角占优的矩阵,这不仅确定了解的唯一性,还保证了迭代法的收敛性。此外, 还可以确定进行 LU 分解,若使用高斯消去法还可以省去选主元的工作。
一、题目重述
解微分方程:
(e y ux ( x, y )) x (e xu y ( x, y )) y ( x y )ux ( x, y ) ( x y )u y ( x, y ) u ( x, y ) y 2e y x 2e x e xy y 2 x 2 1 e xy
0 ⋱ ⋱
0 ⋱ −������������ ,������−1 −������������ ,������−2 ������������ ,������−1
������������1 ������������1 ������������2 ������������2 , ������������ = , ������������ = ⋮ ⋮ ������������ ,������−1 ������������ ,������−1
2.4 模型建立 这里的模型建立就是把差分方程组改写成系数矩阵的形式。经过研究,我觉 得写成如下的系数矩阵不仅看起来简单明了,而且在 matlab 编程时比较方便。 系数矩阵为:Pu=f 其中 P 是(n − 1) 阶方阵,具体如下:
������11 ������12 0 ������11 ������12 ⋱ 0 ������21 ������22 0 ⋱ ⋱ ������11 0 ⋱ ������1,������−1 ������1,������−2 ������1,������−1 ������21 ������22 0 ⋱ ������2,������−1 ������21 ������22 ⋱ 0 ������12 ������13 ⋱ ������1,������−1 0 ⋱ ������1,������−1 ������2,������−2 ������2,������−1 ������������−2,1 ������������−2,2 ⋱ 0 ������������−23 ⋱
5
误差 0.000000021749 0.000000075589 0.000000192565 0.000001495550 0.000001910877 0.000001452420 0.000001013029 0.000001898934
举一个例子:当 i=2,j=3 时,������������������ = ������23 ;当 i=3,j=3 时,������������−1,������ = ������23 。但 是,显然这两个������23 不是同一个数,其大小也不相等。
2.3 差分方程的重新定义 因此,为了避免 2.2 中赋值重复而产生的错误,我在利用 matlab 编程时, 对这些系数变量进行了重新定义:������������������ = ������������������ ,������������������ = ������������ ,������ +1 ,������������������ = ������������ ,������ −1 ,������������������ = ������������ +1,������ ,������������������ = ������������−1,������ .
−2
������������������������ 2 ������������������������ 2 ������������������������ 2 ������������������������ 2
������������ +1/2,������ + ������������−1/2,������ + ������������ ,������ −1/2 + ������������ ,������ +1/2 + ������������������
第二步:
3
bcd1 BCD = 0 第三步: bcd2 ⋱
0 bcdi
������1 ������2 ������2 ������3 ,G = ⋮ ,K = ⋮ ������������−2 ������������−1
P=BCD-diag(G,99)-diag(K,99). P 和 f 的具体构造见附录 6.1 主代码
������������−1,1 ������������−1,2 ⋱ 0
0 ⋱ ⋱
⋱ ������������−1,������−1
������������−1,������−2 ������������−1,������−1
而 u 是(n − 1) 维的列向量,具体如下:
2
2
������ =
������11 ������12 ⋮ ������1,������−1 ������21 ⋮ ⋮ ������������ −1,������−1
4
针对本题,我觉得比较理想的解决方法是采用稀疏矩阵的方式定义 P。 这样可以有效的减小 P 的内存消耗。但是考虑到老师的这次期中作业主要 是考察我们对二阶椭圆偏微分方程的理解与实例操作, 而不是旨在考察我们 的 matlab 编程能力。因此我在此,略作偷懒,把 n 改成 10 或 70(75 以上 内存就不够用了) ,适当的降低精度来得到结果。
已知边界: u(0, y) = 1, u(1, y) = e y , u( x,0) = 1, u( x,1) = e x 求数值解, 把区域 G = [0,1]? [0,1] 分成 h1 = 1/100, h2 = 1/100 ,n=100 注:老师你给的题 F 好像写错了,应该把 y 2ex x2e y 改成 y 2e y x2ex 。
1
������������−1,������ = ������−2 ������������−1/2,������ + ������������ ,������ −1 = ������−2 ������������ ,������ −1/2 + ������������ +1,������ = ������−2 ������������ +1/2,������ − ������������ ,������ +1 = ������−2 ������������ ,������ +1/2 − ������������������ = ������
数值解 1.010050145335 1.020201264438 1.030454341388 1.419066053043 1.491822786765 1.568310733070 1.061837559575 1.127498750514
真实值 1.010050167084 1.020201340027 1.030454533954 1.419067548593 1.49182469764百度文库 1.568312185490 1.061836546545 1.127496851579
四、计算结果
4.1 当 n=10,h=1/10 时的结果 取 n=10,h=1/10 时,matlab 运行的部分结果如图 4.1。表 4.2 为 LU 分解法和高斯 消去法的部分结果 (这两个直接法结果完全一样) , 表 4.3 为迭代法的部分结果。
图 4.1
i,j 1,1 1,2 1,3 5,7 5,8 5,9 6,1 6,2
2
������2,3
0 ⋱ ⋱
0
������������−1,1 ������������−1,2 ������������−1,3 ⋱ ������������−1,������−1
������������−2,,������−1 0
������������−1,1 ������������−1,2 0
图 3.2 3.3.3 问题解决 经过上网查找资料后,我找到了如下几个解决方法。 1)尽量早的分配大的矩阵变量 2)尽量避免产生大的瞬时变量,当它们不用的时候应该及时 clear。 3)尽量的重复使用变量(跟不用的 clear 掉一个意思) 4)将矩阵转化成稀疏形式 5)使用 pack 命令 6)如果可行的话,将一个大的矩阵划分为几个小的矩阵,这样每一次使 用的内存减少。 7)增大内存
《微分方程数值解法》期中作业实验报告
二阶椭圆偏微分方程第一边值问题
姓名: 学号: 班级:
2013 年 11 月 19 日
二阶椭圆偏微分方程第一边值问题 摘要
对于解二阶椭圆偏微分方程第一边值问题, 课本上已经给出了相应的差分方 程。 而留给我的难题就是把差分方程组表示成系数矩阵的形式,以及对系数进行 赋值。解决完这个问题之后,我在利用 matlab 解线性方程组时,又出现“out of memory” 的问题。 因为 99*99 阶的矩阵太大, 超出了分配给 matlab 的使用内存。 退而求其次, 当 n=10, h=1/10 或 n=70, h=1/70 时, 我都得出了很好的计算结果。 然而在解线性方程组时,无论是 LU 分解法或高斯消去法,还是 gauseidel 迭代 法,都能达到很高的精度。 关键字:二阶椭圆偏微分方程差分方程 out of memory LU 分解高斯消去法 gauseidel 迭代法
相关文档
最新文档