2022年二阶椭圆偏微分方程实例求解

1、学习资料收集于网络,仅供参考微分方程数值解法期中作业试验报告二阶椭圆偏微分方程第一边值问题姓名：学号：班级：2022 年 11 月 19 日学习资料学习资料收集于网络,仅供参考二阶椭圆偏微分方程第一边值问题摘要对于解二阶椭圆偏微分方程第一边值问题,课本上已经给出了相应的差分方程；而留给我的难题就是把差分方程组表示成系数矩阵的形式,以及对系数进行 赋值；解决完这个问题之后,我在利用 matlab 解线性方程组时,又显现“out of memory” 的问题；由于 99*99 阶的矩阵太大, 超出了安排给 matlab 的使用内存；退而求其次,当 n=10,h=1/10 或 n=70,h=1/70

2、 时,我都得出了很好的运算结果；然而在解线性方程组时,无论是LU分解法或高斯消去法,仍是gauseidel 迭代法,都能达到很高的精度；关键字：二阶椭圆偏微分方程差分方程out of memory LU分解高斯消去法gauseidel迭代法一、题目重述解微分方程：y e u x , xx e uy , yxy ux , xy uy , u x y , 2 y ey2 xx exy e2 yx21xy e=ex已知边界 :u0,y=1, 1, =ey, ,0=1, ,1求数值解 , 把区域G =0,1.0,1分成h 1=1/100,h2=1/100,n=100注： 老师你给的题 F 似乎写错了,

3、应当把2 y ex2 yx e 改成2 y ey2 xx e ；二、问题分析与模型建立2.1 微分方程上的符号说明2.2 课本上差分方程的缺陷 课本上的差分方程为：学习资料学习资料收集于网络,仅供参考举一个例子：当 i=2,j=3时,；当 i=3,j=3时,；但是,明显这两个 不是同一个数,其大小也不相等；2.3 差分方程的重新定义因此,为了防止 2.2 中赋值重复而产生的错误,我在利用 matlab 编程时,对这些系数变量进行了重新定义：2.4 模型建立这里的模型建立就是把差分方程组改写成系数矩阵的形式；经过讨论,我觉得写成如下的系数矩阵不仅看起来简洁明白,而且在 系数矩阵为： Pu=f其中

4、 P是 阶方阵,详细如下：学习资料matlab 编程时比较便利；学习资料收集于网络,仅供参考而 u 是 维的列向量,详细如下：而 f 是 维的列向量,详细如下：三、求解过程 3.1 对系数矩阵的分析对上述模型的求解就是对线性方程组的求解；通过观看,我发觉 P 是一个对角占优的矩阵,这不仅确定明白的唯独性,仍保证了迭代法的收敛性；此外,仍可以确定进行 LU分解,如使用高斯消去法仍可以省去选主元的工作；3.2matlab 编程由于矩阵阶数过大, 所以此题的编程难点为构造系数矩阵,即对线性方程组的赋值；我采纳的方法是分块赋值；对于 第一步：学习资料P 的赋值,过程如下：学习资料收集于网络,仅供参考其

5、次步：第三步：P=BCD-diagG,99-diagK,99. P 和 f 的详细构造见 附录 6.1 主代码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 左右的内存

6、, 加上其他数据,内存消耗至少在 1G 以上；但是我电脑上安排给 即使在关闭杀毒软件等大部分应用程序后,安排给图 3.2 matlab 的内存只有： 454 MB,matlab 的内存也刚够 1G；3.3.3 问题解决 经过上网查找资料后,我找到了如下几个解决方法；1）尽量早的安排大的矩阵变量 2）尽量防止产生大的瞬时变量,当它们不用的时候应当准时 clear ；3）尽量的重复使用变量（跟不用的 4）将矩阵转化成稀疏形式 5）使用 pack 命令clear 掉一个意思）6）假如可行的话, 将一个大的矩阵划分为几个小的矩阵,这样每一次使学习资料学习资料收集于网络,仅供参考用的内存削减；7）增大内

7、存 针对此题,我觉得比较抱负的解决方法是采纳稀疏矩阵的方式定义 P；这样可以有效的减小P 的内存消耗；但是考虑到老师的这次期中作业主要是考察我们对二阶椭圆偏微分方程的懂得与实例操作,而不是旨在考察我们的 matlab 编程才能；因此我在此,略作偷懒,把n 改成 10 或 70（75 以上内存就不够用了）,适当的降低精度来得到结果；四、运算结果4.1 当 n=10,h=1/10 时的结果 取 n=10,h=1/10 时,matlab 运行的部分结果如图消去法的部分结果（这两个直接法结果完全一样）图 4.1 4.1；表 4.2 为 LU分解法和高斯,表 4.3 为迭代法的部分结果；i,j 数值解真

8、实值误差1,1 1.010050145335 1.010050167084 0.000000021749 1,2 1.020222264438 1.020222240027 0.000000075589 1,3 1.030454341388 1.030454533954 0.000000192565 5,7 1.419066053043 1.419067548593 0.000001495550 5,8 1.491822786765 1.491824697641 0.000001910877 5,9 1.568310733070 1.568312185490 0.000001452420 学习

9、资料学习资料收集于网络,仅供参考6,1 1.061837559575 1.061836546545 0.000001013029 6,2 1.127498750514 1.127496851579 0.000001898934 6,3 1.197219794676 1.197217363122 0.000002431554 6,4 1.271251608972 1.271249150321 0.000002458650 9,7 1.877615353384 1.877610579264 0.000004774120 9,8 2.054437020906 2.054433210644 0.000

10、003810262 9,9 2.247910518362 2.247907986676 0.000002531686 表 4.2 i,j 数值解真实值误差1,1 1.010049929223 1.010050167084 0.000000237861 1,2 1.020222273723 1.020222240027 0.000000466304 1,3 1.030453831277 1.030454533954 0.000000702677 5,5 1.284024086148 1.284025416688 0.000001330540 5,6 1.349856719969 1.349858

11、807576 0.000002087607 5,7 1.419064868769 1.419067548593 0.000002679825 5,8 1.491821975367 1.491824697641 0.000002722274 5,9 1.568310332118 1.568312185490 0.000001853372 6,1 1.061837000239 1.061836546545 0.000000453693 6,2 1.127497727757 1.127496851579 0.000000876177 6,3 1.197218445256 1.197217363122

12、 0.000001082134 9,7 1.877615007879 1.877610579264 0.000004428615 9,8 2.054436783189 2.054433210644 0.000003572545 9,9 2.247910400175 2.247907986676 0.000002413498 表 4.3 4.2 当 n=70,h=1/70 时的结果取 n=70,h=1/70 时,matlab 运行的部分结果如图 4.4（LU分解法）；运算时间为 17多分钟（ 1043 秒）,误差至少在小数点后 9 位；学习资料学习资料收集于网络,仅供参考图 4.4 五、结论分析

13、由于本人的电脑比较破旧, 内存不是很大, 再加上没有实行稀疏矩阵的储备方式,难以达到 n=100,h=1/100 的运算要求；但改为 n=10,h=1/10 或 n=70,h=1/70后,运算精度也非常抱负；特别是后者,其误差至少在小数点后 9 位,在比较使用哪种方法解线性方程组时,直接法中的 LU分解法和高斯消去法的运算结果是完全相等的； 而 gauseidel 迭代法运算结果按个人设定的运算精度而定；对于这种大型线性方程组来说, 迭代法从运算速度和运算机储备方面来看具有超过直接法的打算性优点；对于稀疏方程组来说, 迭代法非常有效； 从此题的实际情形来看,当 n=70 时, LU 分解的直接

14、法的运算时间为 17 分钟左右,而gauseidel 迭代法的运算时间为 20 秒左右,这与以上分析的情形一样； 但是当 n越来越大时 （从 n=10 起）,固定迭代步数的迭代法解的精度越来越差,为了达到与直接法相同的精度, 必需提高迭代步数； 然而这又会加大运算量, 使运算速度变慢（见表 5.1）；所以迭代法是在运算精度和运算速度之间的权衡取舍；n=50,h=1/50 迭代精度迭代步数迭代时间1.00E-04 1555 13.3 秒1.00E-06 2681 21.6 秒1.00E-08 3808 29.8 秒表 5.1 仅从此题运算结果来看（ n=10 时）：当 x,y 靠近 0 时,直接

15、法的数值解更精确,而当 x,y 靠近 1 时,迭代法的数值解更精确；这与 gauseidel 迭代法的算法特点相符合； 由于 gauseidel 迭代法后面的解在迭代时要用到前面的解,所以 x,y 靠近 1 的后面的解更为精确；六、附录6.1 主代码 主代码中解决了对系数矩阵的赋值,即写出了线性方程组； 在解线性方程组时可以选用 LU分解法、高斯消去法和迭代法中的任意一种；function n=Middle_Term_Bymyselft%t为区间划分数 tic; %定义函数及初始化基本变量 FA=x,yexpy; FB=x,yexpx; FC=x,yx+y; FD=x,yx-y; FE=x,y

16、1; FF=x,yy2+x2+1*expx*y-y2*expy+x2*expx*expx*y; FU=x,yexpx*y;%真实值 n=t;%区间划分为 n*n 学习资料学习资料收集于网络,仅供参考h=1/n;%h 为单位区间长度 A=zeros2*n-1; B=zeros2*n-1; C=zerosn-1; D=zerosn-1; E=zerosn-1; F=zerosn-1; U=zerosn-1;%真实值的矩阵表示 u=zerosn-1*n-1,1;%真实值的数列表示 error=zerosn-1*n-1,1;%误差 P=zerosn-1*n-1;%最终要解的方程组的系数矩阵 ,即平常的

17、 A f=zerosn-1*n-1,1; %即平常的 b ff=zerosn-1; %ffi,j=fi-1）*9+j BCD=;%记录三对角部分 G=;%记录上三角的那一斜条 K=;%记录下三角的那一斜条 b=zerosn-1;%表示 ai,j c=zerosn-1;%表示 ai,j+1 d=zerosn-1;%表示 ai,j-1 g=zerosn-1;%表示 ai+1,j k=zerosn-1;%表示 ai-1,j %对函数进行赋值 x=zeros2*n-1,1; y=zeros2*n-1,1; for i=1:2*n-1 %使 A.B下标 i-1/2 变为 2i-1 for j=1:2*n

18、-1 xi=i*h/2; yj=j*h/2; Ai,j=FAxi,yj; Bi,j=FBxi,yj; end end x=zerosn-1,1; y=zerosn-1,1; for i=1:n-1 for j=1:n-1 xi=i*h; yj=j*h; Ci,j=FCxi,yj; Di,j=FDxi,yj; Ei,j=1; Fi,j=FFxi,yj; Ui,j=FUxi,yj; end 学习资料学习资料收集于网络,仅供参考end %对最终要解的方程组的系数矩阵进行赋值 for i=1:n-1 bcd=zerosn-1; bb=; cc=; dd=; gg=; kk=; for j=1:n-1

19、bi,j=A2*i+1,2*j+A2*i-1,2*j+B2*i,2*j+1+B2*i,2*j-1/h2+Ei,j; ci,j=B2*i,2*j+1-h*Di,j/2/h2; di,j=B2*i,2*j-1+h*Di,j/2/h2; gi,j=A2*i+1,2*j-h*Ci,j/2/h2; ki,j=A2*i-1,2*j+h*Ci,j/2/h2; bb=bb bi,j; if j=2 dd=dd di,j; end gg=gg gi,j; kk=kk ki,j; %给 f 赋值 if i=1 ffi,j=Fi,j+ki,j*1;%边值为 1 elseif i=n-1 ffi,j=Fi,j+gi,

20、j*Ai,2*j;%A 中 i 取值无所谓,不影响 else end ffi,j=Fi,j; if j=1 ffi,j=ffi,j+di,j*1;% 边值为 1 elseif j=n-1 ffi,j=ffi,j+ci,j*B2*i,j; end end fi-1*n-1+j,1 = ffi,j;%你应当懂的,坐标变换bcd=diagbb-diagcc,1-diagdd,-1; BCD=blkdiagBCD,bcd; if i=2 K=K kk; end end P=BCD-diagG,n-1-diagK,-n+1; %BCD=BCD-diagG,n-1-diagK,-n+1;x=Doolitt

21、leBCD,f; 这样是不是可以削减点内存 消耗x=DoolittleP,f;%LU分解解方程组,这里也可以输入其他解方程组的方法 %x=GaussXQByOrderP,f % 高斯消去法 %x0=onesn-12,1;x=gauseidelP,f,x0 %gauseidel 迭代法 for i=1:n-1 for j=1:n-1 ui-1*n-1+j=Ui,j; end end error=absx-u; result=x u error; time = toc; disp运算时间为： ; disptime; disp-; format long; disp运算结果为： ; disp 数值解真实值误差 ; dispresult; 6.2LU分解解线性方程组 function x,L,U= Doolittle A,b N = sizeA; n = N1; L = eyen,n; %L的对角元素为 1 U = zerosn,n; U1,1:n = A1,1:n; %U的第一行 L1:n,1 = A1:n,1/U1,1; %L的第一列for k=2:n for i=k:n Uk,i = Ak,i-Lk,1:k-1*U1:k-1,i; %U的第 k 行 end for j=k+1:n 学习资料学习资料收集于网络,仅供参考Lj,k =