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

微分方程數(shù)值解法期中作業(yè)實驗報告二階橢圓偏微分方程第一邊值問題姓名： 學(xué)號：班級：2013年11月19日二階橢圓偏微分方程第一邊值問題摘要對于解二階橢圓偏微分方程第一邊值問題，課本上已經(jīng)給出了相應(yīng)的差分方程。而留給我的難題就是把差分方程組表示成系數(shù)矩陣的形式，以及對系數(shù)進(jìn)行賦值。解決完這個問題之后，我在利用matlab解線性方程組時，又出現(xiàn)“out of memory”的問題。因為99*99階的矩陣太大，超出了分配給matlab的使用內(nèi)存。退而求其次，當(dāng)n=10，h=1/10或n=70，h=1/70時，我都得出了很好的計算結(jié)果。然而在解線性方程組時，無論是LU分解法或高斯消去法，還是gauseidel迭代法，都能達(dá)到很高的精度。 關(guān)鍵字：二階橢圓偏微分方程 差分方程 out of memory LU分解 高斯消去法 gauseidel迭代法一、題目重述解微分方程：已知邊界:求數(shù)值解, 把區(qū)域分成，n=100注：老師你給的題F好像寫錯了，應(yīng)該把改成。二、問題分析與模型建立2.1微分方程上的符號說明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/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舉一個例子：當(dāng)i=2,j=3時，aij=a23；當(dāng)i=3,j=3時，ai-1，j=a23。但是，顯然這兩個a23不是同一個數(shù)，其大小也不相等。2.3差分方程的重新定義因此，為了避免2.2中賦值重復(fù)而產(chǎn)生的錯誤，我在利用matlab編程時，對這些系數(shù)變量進(jìn)行了重新定義：bij=aij，cij=ai,j+1，dij=ai,j-1，gij=ai+1,j，kij=ai-1,j.2.4模型建立 這里的模型建立就是把差分方程組改寫成系數(shù)矩陣的形式。經(jīng)過研究，我覺得寫成如下的系數(shù)矩陣不僅看起來簡單明了，而且在matlab編程時比較方便。系數(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,2000cn-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對系數(shù)矩陣的分析對上述模型的求解就是對線性方程組的求解。通過觀察，我發(fā)現(xiàn)P是一個對角占優(yōu)的矩陣，這不僅確定了解的唯一性，還保證了迭代法的收斂性。此外，還可以確定進(jìn)行LU分解，若使用高斯消去法還可以省去選主元的工作。3.2matlab編程因為矩陣階數(shù)過大，所以此題的編程難點(diǎn)為構(gòu)造系數(shù)矩陣，即對線性方程組的賦值。我采用的方法是分塊賦值。對于P的賦值，過程如下：第一步：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時，運(yùn)行編好的matlab程序時，會出現(xiàn)如圖3.1的錯誤提示。圖3.13.3.2問題分析在matlab的命令窗口輸入“memory”，出現(xiàn)如圖3.2的內(nèi)存使用情況，可以得出：Memory used by MATLAB: 454 MB (4.760e+008 bytes)。，若不用稀疏矩陣定義P，經(jīng)過粗略計算，我發(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)查找資料后，我找到了如下幾個解決方法。1）盡量早的分配大的矩陣變量2）盡量避免產(chǎn)生大的瞬時變量，當(dāng)它們不用的時候應(yīng)該及時clear。3）盡量的重復(fù)使用變量（跟不用的clear掉一個意思）4）將矩陣轉(zhuǎn)化成稀疏形式5）使用pack命令6）如果可行的話，將一個大的矩陣劃分為幾個小的矩陣，這樣每一次使用的內(nèi)存減少。7）增大內(nèi)存針對本題，我覺得比較理想的解決方法是采用稀疏矩陣的方式定義P。這樣可以有效的減小P的內(nèi)存消耗。但是考慮到老師的這次期中作業(yè)主要是考察我們對二階橢圓偏微分方程的理解與實例操作，而不是旨在考察我們的matlab編程能力。因此我在此，略作偷懶，把n改成10或70（75以上內(nèi)存就不夠用了），適當(dāng)?shù)慕档途葋淼玫浇Y(jié)果。四、計算結(jié)果4.1當(dāng) n=10,h=1/10時的結(jié)果取n=10,h=1/10時，matlab運(yùn)行的部分結(jié)果如圖4.1。表4.2為LU分解法和高斯消去法的部分結(jié)果（這兩個直接法結(jié)果完全一樣），表4.3為迭代法的部分結(jié)果。圖4.1i,j數(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.568312185490 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.000003810262 9,92.247910518362 2.247907986676 0.000002531686 表4.2i,j數(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 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.000001082134 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時的結(jié)果取n=70,h=1/70時，matlab運(yùn)行的部分結(jié)果如圖4.4（LU分解法）。計算時間為17多分鐘（1043秒），誤差至少在小數(shù)點(diǎn)后9位。圖4.4五、結(jié)論分析由于本人的電腦比較破舊，內(nèi)存不是很大，再加上沒有采取稀疏矩陣的存儲方式，難以達(dá)到n=100，h=1/100的計算要求。但改為n=10，h=1/10或n=70，h=1/70后，計算精度也十分理想。尤其是后者，其誤差至少在小數(shù)點(diǎn)后9位， 在比較使用哪種方法解線性方程組時，直接法中的LU分解法和高斯消去法的計算結(jié)果是完全相等的。而gauseidel迭代法計算結(jié)果按個人設(shè)定的計算精度而定。對于這種大型線性方程組來說，迭代法從計算速度和計算機(jī)存儲方面來看具有超過直接法的決定性優(yōu)點(diǎn)。對于稀疏方程組來說，迭代法十分有效。從本題的實際情況來看，當(dāng)n=70時，LU分解的直接法的計算時間為17分鐘左右，而gauseidel迭代法的計算時間為20秒左右，這與以上分析的情況一致。但是當(dāng)n越來越大時（從n=10起），固定迭代步數(shù)的迭代法解的精度越來越差，為了達(dá)到與直接法相同的精度，必須提高迭代步數(shù)。然而這又會加大計算量，使計算速度變慢（見表5.1）。所以迭代法是在計算精度和計算速度之間的權(quán)衡取舍。n=50，h=1/50迭代精度迭代步數(shù)迭代時間1.00E-04155513.3秒1.00E-06268121.6秒1.00E-08380829.8秒表5.1僅從本題計算結(jié)果來看（n=10時）：當(dāng)x，y靠近0時，直接法的數(shù)值解更準(zhǔn)確，而當(dāng)x，y靠近1時，迭代法的數(shù)值解更準(zhǔn)確。這與gauseidel迭代法的算法特點(diǎn)相符合。因為gauseidel迭代法后面的解在迭代時要用到前面的解，所以x，y靠近1的后面的解更為精確。六、附錄6.1主代碼主代碼中解決了對系數(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);%真實值n=t;%區(qū)間劃分為n*nh=1/n;%h為單位區(qū)間長度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);%真實值的矩陣表示u=zeros(n-1)*(n-1),1);%真實值的數(shù)列表示error=zeros(n-1)*(n-1),1);%誤差P=zeros(n-1)*(n-1);%最終要解的方程組的系數(shù)矩陣,即平時的Af=zeros(n-1)*(n-1),1); %即平時的bff=zeros(n-1); %ff(i,j)=f(i-1）*9+j)BCD=;%記錄三對角部分G=;%記錄上三角的那一斜條K=;%記錄下三角的那一斜條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)%對函數(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,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%對最終要解的方程組的系數(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*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=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)*1;%邊值為1 elseif i=n-1 ff(i,j)=F(i,j)+g(i,j)*A(i,2*j);%A中i取值無所謂，不影響 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=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,j); endenderror=abs(x-u);result=x u error;time = toc;disp(計算時間為：);disp(time);disp(-);format long;disp(計算結(jié)果為：);disp( 數(shù)值解 真實值 誤差);disp(result); 6.2LU分解解線性方程組function x,L,U= Doolittle (A,b)N = size(A);n = N(1);L = eye(n,n); %L的對角元素為1U = zeros(n,n);U(1,1:n) = A(1,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) = (