二階橢圓偏微分方程實(shí)例求解(附matlab代碼)

1、微分方程數(shù)值解法期中作業(yè)實(shí)驗(yàn)報(bào)告二階橢圓偏微分方程第一邊值問題姓名： 學(xué)號(hào)：班級(jí)：2013年11月19日二階橢圓偏微分方程第一邊值問題摘要對(duì)于解二階橢圓偏微分方程第一邊值問題，課本上已經(jīng)給出了相應(yīng)的差分方程。而留給我的難題就是把差分方程組表示成系數(shù)矩陣的形式，以及對(duì)系數(shù)進(jìn)行賦值。解決完這個(gè)問題之后，我在利用matlab解線性方程組時(shí)，又出現(xiàn)“out of memory”的問題。因?yàn)?9*99階的矩陣太大，超出了分配給matlab的使用內(nèi)存。退而求其次，當(dāng)n=10，h=1/10或n=70，h=1/70時(shí)，我都得出了很好的計(jì)算結(jié)果。然而在解線性方程組時(shí)，無(wú)論是LU分解法或高斯消去法，還是gause

2、idel迭代法，都能達(dá)到很高的精度。 關(guān)鍵字：二階橢圓偏微分方程 差分方程 out of memory LU分解 高斯消去法 gauseidel迭代法一、題目重述解微分方程：已知邊界:求數(shù)值解, 把區(qū)域分成，n=100注：老師你給的題F好像寫錯(cuò)了，應(yīng)該把改成。二、問題分析與模型建立2.1微分方程上的符號(hào)說明Ax,y=ey Bx,y=ex Cx,y=x+y Dx,y=x-y Ex,y=1 Fx,y=2.2課本上差分方程的缺陷課本上的差分方程為：aijuij-ai-1,jui-1,j+ai,j-1ui,j-1+ai+1,jui+1,j+ai,j+1ui,j+1=Fijai-1,j=h-2Ai-1/

3、2,j+hCij2ai,j-1=h-2Bi,j-1/2+hDij2ai+1,j=h-2Ai+1/2,j-hCij2ai,j+1=h-2Bi,j+1/2-hDij2aij=h-2Ai+1/2,j+Ai-1/2,j+Bi,j-1/2+Bi,j+1/2+Eij舉一個(gè)例子：當(dāng)i=2,j=3時(shí)，aij=a23；當(dāng)i=3,j=3時(shí)，ai-1，j=a23。但是，顯然這兩個(gè)a23不是同一個(gè)數(shù)，其大小也不相等。2.3差分方程的重新定義因此，為了避免2.2中賦值重復(fù)而產(chǎn)生的錯(cuò)誤，我在利用matlab編程時(shí)，對(duì)這些系數(shù)變量進(jìn)行了重新定義：bij=aij，cij=ai,j+1，dij=ai,j-1，gij=ai+1,

4、j，kij=ai-1,j.2.4模型建立 這里的模型建立就是把差分方程組改寫成系數(shù)矩陣的形式。經(jīng)過研究，我覺得寫成如下的系數(shù)矩陣不僅看起來簡(jiǎn)單明了，而且在matlab編程時(shí)比較方便。系數(shù)矩陣為：Pu=f其中P是（n-1）2階方陣，具體如下：b11c110d12b12000c1,n-2d1,n-1b1,n-1g11g12g13g1,n-1k21k22k2,3k2,n-1b21c210d22b22000c2,n-2d1,n-1b2,n-100gn-2,1gn-2,2gn-23gn-2,n-1kn-1,1kn-1,2kn-1,3kn-1,n-1bn-1,1cn-1,10dn-1,2bn-1,2000

5、cn-1,n-2dn-1,n-1bn-1,n-1而u是（n-1）2維的列向量，具體如下：u=u11u12u1,n-1u21un-1,n-1而f是（n-1）2維的列向量，具體如下：f=f11f12f1,n-1f21fn-1,n-1三、求解過程3.1對(duì)系數(shù)矩陣的分析對(duì)上述模型的求解就是對(duì)線性方程組的求解。通過觀察，我發(fā)現(xiàn)P是一個(gè)對(duì)角占優(yōu)的矩陣，這不僅確定了解的唯一性，還保證了迭代法的收斂性。此外，還可以確定進(jìn)行LU分解，若使用高斯消去法還可以省去選主元的工作。3.2matlab編程因?yàn)榫仃囯A數(shù)過大，所以此題的編程難點(diǎn)為構(gòu)造系數(shù)矩陣，即對(duì)線性方程組的賦值。我采用的方法是分塊賦值。對(duì)于P的賦值，過程如

6、下：第一步：bcdi=bi1-ci10-di2bi2000-ci,n-2-di,n-1bi,n-1,gi=gi1gi2gi,n-1,ki=ki1ki2ki,n-1第二步：BCD=bcd1bcd200bcdi ,G=g1g2gn-2,K=k2k3kn-1第三步： P=BCD-diag(G,99)-diag(K,99).P和 f的具體構(gòu)造見附錄6.1主代碼3.3編程求解過程中的問題 3.3.1問題產(chǎn)生當(dāng)按照老師要求，n=100，h=1/100時(shí)，運(yùn)行編好的matlab程序時(shí)，會(huì)出現(xiàn)如圖3.1的錯(cuò)誤提示。圖3.13.3.2問題分析在matlab的命令窗口輸入“memory”，出現(xiàn)如圖3.2的內(nèi)存使用

7、情況，可以得出：Memory used by MATLAB: 454 MB (4.760e+008 bytes)。，若不用稀疏矩陣定義P，經(jīng)過粗略計(jì)算，我發(fā)現(xiàn)矩陣P就要占800MB左右的內(nèi)存，加上其他數(shù)據(jù)，內(nèi)存消耗至少在1G以上。但是我電腦上分配給matlab的內(nèi)存只有：454 MB，即使在關(guān)閉殺毒軟件等大部分應(yīng)用程序后，分配給matlab的內(nèi)存也剛夠1G。圖3.2 3.3.3問題解決經(jīng)過上網(wǎng)查找資料后，我找到了如下幾個(gè)解決方法。1）盡量早的分配大的矩陣變量2）盡量避免產(chǎn)生大的瞬時(shí)變量，當(dāng)它們不用的時(shí)候應(yīng)該及時(shí)clear。3）盡量的重復(fù)使用變量（跟不用的clear掉一個(gè)意思）4）將矩陣轉(zhuǎn)化成稀

8、疏形式5）使用pack命令6）如果可行的話，將一個(gè)大的矩陣劃分為幾個(gè)小的矩陣，這樣每一次使用的內(nèi)存減少。7）增大內(nèi)存針對(duì)本題，我覺得比較理想的解決方法是采用稀疏矩陣的方式定義P。這樣可以有效的減小P的內(nèi)存消耗。但是考慮到老師的這次期中作業(yè)主要是考察我們對(duì)二階橢圓偏微分方程的理解與實(shí)例操作，而不是旨在考察我們的matlab編程能力。因此我在此，略作偷懶，把n改成10或70（75以上內(nèi)存就不夠用了），適當(dāng)?shù)慕档途葋淼玫浇Y(jié)果。四、計(jì)算結(jié)果4.1當(dāng) n=10,h=1/10時(shí)的結(jié)果取n=10,h=1/10時(shí)，matlab運(yùn)行的部分結(jié)果如圖4.1。表4.2為L(zhǎng)U分解法和高斯消去法的部分結(jié)果（這兩個(gè)直接法

9、結(jié)果完全一樣），表4.3為迭代法的部分結(jié)果。圖4.1i,j數(shù)值解真實(shí)值誤差1，11.010050145335 1.010050167084 0.000000021749 1，21.020201264438 1.020201340027 0.000000075589 1,31.030454341388 1.030454533954 0.000000192565 5,71.419066053043 1.419067548593 0.000001495550 5,81.491822786765 1.491824697641 0.000001910877 5,91.568310733070 1.568

10、312185490 0.000001452420 6,11.061837559575 1.061836546545 0.000001013029 6,21.127498750514 1.127496851579 0.000001898934 6,31.197219794676 1.197217363122 0.000002431554 6,41.271251608972 1.271249150321 0.000002458650 9,71.877615353384 1.877610579264 0.000004774120 9,82.054437020906 2.054433210644 0.

11、000003810262 9,92.247910518362 2.247907986676 0.000002531686 表4.2i,j數(shù)值解真實(shí)值誤差1，11.010049929223 1.010050167084 0.000000237861 1，21.020200873723 1.020201340027 0.000000466304 1,31.030453831277 1.030454533954 0.000000702677 5,51.284024086148 1.284025416688 0.000001330540 5,61.349856719969 1.349858807576

12、 0.000002087607 5,71.419064868769 1.419067548593 0.000002679825 5,81.491821975367 1.491824697641 0.000002722274 5,91.568310332118 1.568312185490 0.000001853372 6,11.061837000239 1.061836546545 0.000000453693 6,21.127497727757 1.127496851579 0.000000876177 6,31.197218445256 1.197217363122 0.000001082

13、134 9,71.877615007879 1.877610579264 0.000004428615 9,82.054436783189 2.054433210644 0.000003572545 9,92.247910400175 2.247907986676 0.000002413498 表4.34.2當(dāng) n=70,h=1/70時(shí)的結(jié)果取n=70,h=1/70時(shí)，matlab運(yùn)行的部分結(jié)果如圖4.4（LU分解法）。計(jì)算時(shí)間為17多分鐘（1043秒），誤差至少在小數(shù)點(diǎn)后9位。圖4.4五、結(jié)論分析由于本人的電腦比較破舊，內(nèi)存不是很大，再加上沒有采取稀疏矩陣的存儲(chǔ)方式，難以達(dá)到n=100，h=

14、1/100的計(jì)算要求。但改為n=10，h=1/10或n=70，h=1/70后，計(jì)算精度也十分理想。尤其是后者，其誤差至少在小數(shù)點(diǎn)后9位， 在比較使用哪種方法解線性方程組時(shí)，直接法中的LU分解法和高斯消去法的計(jì)算結(jié)果是完全相等的。而gauseidel迭代法計(jì)算結(jié)果按個(gè)人設(shè)定的計(jì)算精度而定。對(duì)于這種大型線性方程組來說，迭代法從計(jì)算速度和計(jì)算機(jī)存儲(chǔ)方面來看具有超過直接法的決定性優(yōu)點(diǎn)。對(duì)于稀疏方程組來說，迭代法十分有效。從本題的實(shí)際情況來看，當(dāng)n=70時(shí)，LU分解的直接法的計(jì)算時(shí)間為17分鐘左右，而gauseidel迭代法的計(jì)算時(shí)間為20秒左右，這與以上分析的情況一致。但是當(dāng)n越來越大時(shí)（從n=10起

15、），固定迭代步數(shù)的迭代法解的精度越來越差，為了達(dá)到與直接法相同的精度，必須提高迭代步數(shù)。然而這又會(huì)加大計(jì)算量，使計(jì)算速度變慢（見表5.1）。所以迭代法是在計(jì)算精度和計(jì)算速度之間的權(quán)衡取舍。n=50，h=1/50迭代精度迭代步數(shù)迭代時(shí)間1.00E-04155513.3秒1.00E-06268121.6秒1.00E-08380829.8秒表5.1僅從本題計(jì)算結(jié)果來看（n=10時(shí)）：當(dāng)x，y靠近0時(shí)，直接法的數(shù)值解更準(zhǔn)確，而當(dāng)x，y靠近1時(shí)，迭代法的數(shù)值解更準(zhǔn)確。這與gauseidel迭代法的算法特點(diǎn)相符合。因?yàn)間auseidel迭代法后面的解在迭代時(shí)要用到前面的解，所以x，y靠近1的后面的解更為精

16、確。六、附錄6.1主代碼主代碼中解決了對(duì)系數(shù)矩陣的賦值，即寫出了線性方程組。在解線性方程組時(shí)可以選用LU分解法、高斯消去法和迭代法中的任意一種。function n=Middle_Term_Bymyself(t)%t為區(qū)間劃分?jǐn)?shù)tic;%定義函數(shù)及初始化基本變量FA=(x,y)exp(y);FB=(x,y)exp(x);FC=(x,y)x+y;FD=(x,y)x-y;FE=(x,y)1;FF=(x,y)(y2+x2+1)*exp(x*y)-(y2*exp(y)+x2*exp(x)*exp(x*y);FU=(x,y)exp(x*y);%真實(shí)值n=t;%區(qū)間劃分為n*nh=1/n;%h為單位區(qū)間長(zhǎng)

17、度A=zeros(2*n-1);B=zeros(2*n-1);C=zeros(n-1);D=zeros(n-1);E=zeros(n-1);F=zeros(n-1);U=zeros(n-1);%真實(shí)值的矩陣表示u=zeros(n-1)*(n-1),1);%真實(shí)值的數(shù)列表示error=zeros(n-1)*(n-1),1);%誤差P=zeros(n-1)*(n-1);%最終要解的方程組的系數(shù)矩陣,即平時(shí)的Af=zeros(n-1)*(n-1),1); %即平時(shí)的bff=zeros(n-1); %ff(i,j)=f(i-1）*9+j)BCD=;%記錄三對(duì)角部分G=;%記錄上三角的那一斜條K=;%記

18、錄下三角的那一斜條b=zeros(n-1);%表示a(i,j)c=zeros(n-1);%表示a(i,j+1)d=zeros(n-1);%表示a(i,j-1)g=zeros(n-1);%表示a(i+1,j)k=zeros(n-1);%表示a(i-1,j)%對(duì)函數(shù)進(jìn)行賦值x=zeros(2*n-1,1);y=zeros(2*n-1,1);for i=1:2*n-1 %使A.B下標(biāo)i-1/2變?yōu)?i-1 for j=1:2*n-1 x(i)=i*h/2; y(j)=j*h/2; A(i,j)=FA(x(i),y(j); B(i,j)=FB(x(i),y(j); endendx=zeros(n-1,

19、1);y=zeros(n-1,1);for i=1:n-1 for j=1:n-1 x(i)=i*h; y(j)=j*h; C(i,j)=FC(x(i),y(j); D(i,j)=FD(x(i),y(j); E(i,j)=1; F(i,j)=FF(x(i),y(j); U(i,j)=FU(x(i),y(j); endend%對(duì)最終要解的方程組的系數(shù)矩陣進(jìn)行賦值for i=1:n-1 bcd=zeros(n-1); bb=; cc=; dd=; gg=; kk=; for j=1:n-1 b(i,j)=(A(2*i+1,2*j)+A(2*i-1,2*j)+B(2*i,2*j+1)+B(2*i,2

20、*j-1)/h2+E(i,j); c(i,j)=(B(2*i,2*j+1)-h*D(i,j)/2)/h2; d(i,j)=(B(2*i,2*j-1)+h*D(i,j)/2)/h2; g(i,j)=(A(2*i+1,2*j)-h*C(i,j)/2)/h2; k(i,j)=(A(2*i-1,2*j)+h*C(i,j)/2)/h2; bb=bb b(i,j); if j<=n-2 cc=cc c(i,j); end if j>=2 dd=dd d(i,j); end gg=gg g(i,j); kk=kk k(i,j); %給f賦值 if i=1 ff(i,j)=F(i,j)+k(i,j

21、)*1;%邊值為1 elseif i=n-1 ff(i,j)=F(i,j)+g(i,j)*A(i,2*j);%A中i取值無(wú)所謂，不影響 else ff(i,j)=F(i,j); end if j=1 ff(i,j)=ff(i,j)+d(i,j)*1;%邊值為1 elseif j=n-1 ff(i,j)=ff(i,j)+c(i,j)*B(2*i,j); end f(i-1)*(n-1)+j,1) = ff(i,j);%你應(yīng)該懂的，坐標(biāo)變換 end bcd=diag(bb)-diag(cc,1)-diag(dd,-1); BCD=blkdiag(BCD,bcd); if i<=n-2 G=G

22、 gg; end if i>=2 K=K kk; end endP=BCD-diag(G,n-1)-diag(K,-n+1);%BCD=BCD-diag(G,n-1)-diag(K,-n+1);x=Doolittle(BCD,f);這樣是不是可以減少點(diǎn)內(nèi)存消耗x=Doolittle(P,f);%LU分解解方程組，這里也可以輸入其他解方程組的方法%x=GaussXQByOrder(P,f) %高斯消去法%x0=ones(n-1)2,1);x=gauseidel(P,f,x0) %gauseidel迭代法for i=1:n-1 for j=1:n-1 u(i-1)*(n-1)+j)=U(i,

23、j); endenderror=abs(x-u);result=x u error;time = toc;disp('計(jì)算時(shí)間為：');disp(time);disp('-');format long;disp('計(jì)算結(jié)果為：');disp(' 數(shù)值解 真實(shí)值 誤差');disp(result); 6.2LU分解解線性方程組function x,L,U= Doolittle (A,b)N = size(A);n = N(1);L = eye(n,n); %L的對(duì)角元素為1U = zeros(n,n);U(1,1:n) = A(1,

24、1:n); %U的第一行L(1:n,1) = A(1:n,1)/U(1,1); %L的第一列for k=2:n for i=k:n U(k,i) = A(k,i)-L(k,1:(k-1)*U(1:(k-1),i); %U的第k行 end for j=(k+1):n L(j,k) = (A(j,k)-L(j,1:(k-1)*U(1:(k-1),k)/U(k,k); %L的第k列 endendy = SolveDownTriangle(L,b);x = SolveUpTriangle(U,y); %求解方程6.3高斯消去法function x,XA=GaussXQByOrder(A,b)N = size(A);n = N(1);for i=1:(n-1) for j=(i+1):n if(A(i,i)=0) disp('對(duì)角元素為0！'); %防止對(duì)角元素為0 return; end l = A(j,i); m = A(i,i); A(j,1:n)=A(j,1:n)