迭代法求微分方程

1、迭代法的应用拟合求解偏微分方程摘要：现实生活中，因为各种各样的因素很难完全考虑在内使得很多问题很难得到精确解，实际上也不需要得到精确解。例如普通工程中可以允许5%的误差。如果我们能够在规定的误差内较为方便的得到近似解，也是一个很好的解决方式。数值分析中介绍的一些迭代法例如Jacobi迭代和Gauss-Seidel迭代可以较好的解决这些问题，本文中将以这两种迭代方式为主解一偏微分方程，而该偏微分方程来自于现实中的洋流问题。关键词：Jacobi迭代 Gauss-Seidel迭代 近似解一、问题重述研究一地球表面的一正方形洋流，其中正方向代表地球的正东，正方向代表地球的正北方向。洋流由一流函数定义，

2、从而实际流速可以由向量。如果地球是平的，流函数应该满足拉普拉斯方程，考虑地球的曲率在内，由纬度决定的科氏力应该考虑在内，得到一对流扩散的方程：考虑到在其边界没有流动，添加边界条件。这样得到的解为在任一点。事实上还应注意到海洋表面风的存在，因此在方程右端添加一非平凡解正弦强迫项。最终得到满足的方程：边界上满足：。要求求解出。二、解决方案理论计算是很繁琐的，将与有关。不妨换一种思路，使用迭代拟合的方法并利用计算机强大的运算能力得到其误差允许范围内的近似解。具体实施如下：将分为格，在正方形区域所有点都成立转化为在所有格点上都成立，进行数值拟合。由，可得将上面三个近似公式带入原方程，得下面分别使用Ja

3、cobi迭代和Gauss-Seidel迭代编得程序，来近似求解该问题，并讨论在这两种迭代方法下，的变化对迭代步数以及迭代结果的影响。（1）Jacobi迭代（程序见附录）当并设置迭代精度时，迭代步数iter=51，绘出得到的近似解的等高线如下所示（本文中所有图，为排版方便原图缩小了一定比例）：下面讨论在Jacobi迭代方法下，的变化对迭代步数以及迭代结果的影响。 固定和精度，增加观察现象n=10时 ,迭代步数iter=80，绘出得到的近似解的等高线如下所示：n=15时 ,迭代步数iter=174，绘出得到的近似解的等高线如下所示： n=30时 ,迭代步数iter=645，绘出得到的近似解的等高线

4、如下所示：n=54时 ,迭代步数iter=1939，绘出得到的近似解的等高线如下所示：n=55时 ,程序运行显示迭代次数超过了限制次数，无法得到解。可见固定和精度，随着的增大，迭代步数不断增加，得到的解也更加趋于真实解。 固定和精度，减小观察现象时 ,迭代步数iter=41，绘出得到的近似解的等高线如下所示：时 ,迭代步数iter=33，绘出得到的近似解的等高线如下所示：时 ,迭代步数iter=1942，不收敛。可见固定和精度，随着的减小，迭代步数不断减少，得到的解的等高线也越来越陡，不再呈一类似圆形。（2）Gauss-Seidel迭代（程序见附录）当并设置迭代精度时，迭代步数iter=31，

5、绘出得到的近似解的等高线如下所示：下面讨论在Gauss-Seidel迭代方法下，的变化对迭代步数以及迭代结果的影响。 固定和精度，增加观察现象n=10时 ,迭代步数iter=46，绘出得到的近似解的等高线如下所示：n=15时 ,迭代步数iter=47，绘出得到的近似解的等高线如下所示：n=30时 ,迭代步数iter=360，绘出得到的近似解的等高线如下所示：n=54 时 ,迭代步数iter=1076，绘出得到的近似解的等高线如下所示： n=55 时 ,迭代步数iter=1113，绘出得到的近似解的等高线如下所示：与Jacobi迭代类似，固定和精度，随着的增大，迭代步数不断增加，得到的解也更加趋

6、于真实解。这里要指出与Jacobi迭代不同的一点是，在相同的输入参数下，GS迭代显然迭代步数少的多，也说明了GS迭代的优越性。（3）最后假设地球旋转反向，即原程序中反号，当并设置迭代精度时，观察其运行结果可见等高线中间圆形偏向于右方，这可以由地球偏转方向相反，科氏力反向解释。三、小结数值计算方法在解决实际问题时是非常有用的，尤其是对于理论计算很难求得且不需要过于精确值时，使用迭代方法可以求得误差允许范围内的近似解。四、附录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