迭代法求微分方程

1、迭代法的應(yīng)用擬合求解偏微分方程摘要：現(xiàn)實(shí)生活中，因?yàn)楦鞣N各樣的因素很難完全考慮在內(nèi)使得很多問題很難得到精確解，實(shí)際上也不需要得到精確解。例如普通工程中可以允許5%的誤差。如果我們能夠在規(guī)定的誤差內(nèi)較為方便的得到近似解，也是一個(gè)很好的解決方式。數(shù)值分析中介紹的一些迭代法例如Jacobi迭代和Gauss-Seidel迭代可以較好的解決這些問題，本文中將以這兩種迭代方式為主解一偏微分方程，而該偏微分方程來自于現(xiàn)實(shí)中的洋流問題。關(guān)鍵詞：Jacobi迭代 Gauss-Seidel迭代 近似解一、問題重述研究一地球表面的一正方形洋流，其中正方向代表地球的正東，正方向代表地球的正北方向。洋流由一流函數(shù)定義，

2、從而實(shí)際流速可以由向量。如果地球是平的，流函數(shù)應(yīng)該滿足拉普拉斯方程，考慮地球的曲率在內(nèi)，由緯度決定的科氏力應(yīng)該考慮在內(nèi)，得到一對(duì)流擴(kuò)散的方程：考慮到在其邊界沒有流動(dòng)，添加邊界條件。這樣得到的解為在任一點(diǎn)。事實(shí)上還應(yīng)注意到海洋表面風(fēng)的存在，因此在方程右端添加一非平凡解正弦強(qiáng)迫項(xiàng)。最終得到滿足的方程：邊界上滿足：。要求求解出。二、解決方案理論計(jì)算是很繁瑣的，將與有關(guān)。不妨換一種思路，使用迭代擬合的方法并利用計(jì)算機(jī)強(qiáng)大的運(yùn)算能力得到其誤差允許范圍內(nèi)的近似解。具體實(shí)施如下：將分為格，在正方形區(qū)域所有點(diǎn)都成立轉(zhuǎn)化為在所有格點(diǎn)上都成立，進(jìn)行數(shù)值擬合。由，可得將上面三個(gè)近似公式帶入原方程，得下面分別使用Ja

3、cobi迭代和Gauss-Seidel迭代編得程序，來近似求解該問題，并討論在這兩種迭代方法下，的變化對(duì)迭代步數(shù)以及迭代結(jié)果的影響。（1）Jacobi迭代（程序見附錄）當(dāng)并設(shè)置迭代精度時(shí)，迭代步數(shù)iter=51，繪出得到的近似解的等高線如下所示（本文中所有圖，為排版方便原圖縮小了一定比例）：下面討論在Jacobi迭代方法下，的變化對(duì)迭代步數(shù)以及迭代結(jié)果的影響。 固定和精度，增加觀察現(xiàn)象n=10時(shí) ,迭代步數(shù)iter=80，繪出得到的近似解的等高線如下所示：n=15時(shí) ,迭代步數(shù)iter=174，繪出得到的近似解的等高線如下所示： n=30時(shí) ,迭代步數(shù)iter=645，繪出得到的近似解的等高線

4、如下所示：n=54時(shí) ,迭代步數(shù)iter=1939，繪出得到的近似解的等高線如下所示：n=55時(shí) ,程序運(yùn)行顯示迭代次數(shù)超過了限制次數(shù)，無法得到解?？梢姽潭ê途龋S著的增大，迭代步數(shù)不斷增加，得到的解也更加趨于真實(shí)解。 固定和精度，減小觀察現(xiàn)象時(shí) ,迭代步數(shù)iter=41，繪出得到的近似解的等高線如下所示：時(shí) ,迭代步數(shù)iter=33，繪出得到的近似解的等高線如下所示：時(shí) ,迭代步數(shù)iter=1942，不收斂?？梢姽潭ê途?，隨著的減小，迭代步數(shù)不斷減少，得到的解的等高線也越來越陡，不再呈一類似圓形。（2）Gauss-Seidel迭代（程序見附錄）當(dāng)并設(shè)置迭代精度時(shí)，迭代步數(shù)iter=31，

5、繪出得到的近似解的等高線如下所示：下面討論在Gauss-Seidel迭代方法下，的變化對(duì)迭代步數(shù)以及迭代結(jié)果的影響。 固定和精度，增加觀察現(xiàn)象n=10時(shí) ,迭代步數(shù)iter=46，繪出得到的近似解的等高線如下所示：n=15時(shí) ,迭代步數(shù)iter=47，繪出得到的近似解的等高線如下所示：n=30時(shí) ,迭代步數(shù)iter=360，繪出得到的近似解的等高線如下所示：n=54 時(shí) ,迭代步數(shù)iter=1076，繪出得到的近似解的等高線如下所示： n=55 時(shí) ,迭代步數(shù)iter=1113，繪出得到的近似解的等高線如下所示：與Jacobi迭代類似，固定和精度，隨著的增大，迭代步數(shù)不斷增加，得到的解也更加趨

6、于真實(shí)解。這里要指出與Jacobi迭代不同的一點(diǎn)是，在相同的輸入?yún)?shù)下，GS迭代顯然迭代步數(shù)少的多，也說明了GS迭代的優(yōu)越性。（3）最后假設(shè)地球旋轉(zhuǎn)反向，即原程序中反號(hào)，當(dāng)并設(shè)置迭代精度時(shí)，觀察其運(yùn)行結(jié)果可見等高線中間圓形偏向于右方，這可以由地球偏轉(zhuǎn)方向相反，科氏力反向解釋。三、小結(jié)數(shù)值計(jì)算方法在解決實(shí)際問題時(shí)是非常有用的，尤其是對(duì)于理論計(jì)算很難求得且不需要過于精確值時(shí)，使用迭代方法可以求得誤差允許范圍內(nèi)的近似解。四、附錄1.Jacobi迭代function phi,res=fast_current(n,epsilon,tol)% Jacobi iteration for ocean curr

7、ent model.%phi,res = fast_current(n,epsilon,tol);%input parameters:% n : subdivision parameter%epsilon : coriolis force parameter% tol : iteration tolerance%output parameters:% phi : solution (matrix of mesh point values)% res : iteration convergence recordh=1/n;fprintf(subdivision parameter : %gn,h

8、)fprintf(cell peclet number : %gn,h/(2*epsilon)x=0:h:1; y=x;phi=zeros(n+1,n+1);new_phi=zeros(n+1,n+1);% set up rhs vectortau=ones(n+1,1)*sin(h*pi*0:n);it=0;diff=inf;% loop while iteration not convergedticwhile difftol it=it+1; if it1999, error(maximum number of iterations exceeded!) end diff=0; mult

9、=h*h/(4*epsilon); re=h/(2*epsilon); for i=2:n for j=2:n new_phi(i,j)=(phi(i+1,j)+phi(i-1,j)+phi(i,j+1)+phi(i,j-1)/4 . +re*(phi(i+1,j)-phi(i-1,j)/4)+mult*tau(i,j); end end diff=norm(new_phi-phi,inf); phi=new_phi; res(it)=diff;endtt=toc;fprintf(iteration converged in %g iterationsn,it)fprintf(and took

10、 %g secondsn,tt)% plot solution and convergence historyfigure(1)semilogy(1:it,res,b-)axis(square)title(Jacobi convergence history)figure(2)contour(phi,18)axis(ij),axis(off),axis(square)title(integrated current : n is ,num2str(n), . :epsilon is ,num2str(epsilon)return2.GS迭代function phi,res=fast_curre

11、ntgs(n,epsilon,tol)% GS iteration for ocean current model.%phi,res = fast_current(n,epsilon,tol);%input parameters:% n : subdivision parameter%epsilon : coriolis force parameter% tol : iteration tolerance%output parameters:% phi : solution (matrix of mesh point values)% res : iteration convergence r

12、ecordh=1/n;fprintf(subdivision parameter : %gn,h)fprintf(cell peclet number : %gn,h/(2*epsilon)x=0:h:1; y=x;phi=zeros(n+1,n+1);new_phi=zeros(n+1,n+1);% set up rhs vectortau=ones(n+1,1)*sin(h*pi*0:n);it=0;diff=inf;% loop while iteration not convergedticwhile difftol it=it+1; if it1999, error(maximum

13、number of iterations exceeded!) end diff=0; mult=h*h/(4*epsilon); re=h/(2*epsilon); old_phi=phi; for i=2:n for j=2:n phi(i,j)=(phi(i+1,j)+phi(i-1,j)+phi(i,j+1)+phi(i,j-1)/4 . +re*(phi(i+1,j)-phi(i-1,j)/4)+mult*tau(i,j); end end diff=norm(phi-old_phi,inf); %phi=new_phi; res(it)=diff;endtt=toc;fprintf(iteration converged in %g iterationsn,it)fprintf(and took %g secondsn,tt)% plot solution and convergence hist