数理方程基于matlab的数值解法

数理方程数值解法与其在matlab软件上的实现张体强1026222廖荣发1026226[摘要]数学物理方程的数值解在实际生活中越来越使用，首先基于偏微分数值解的思想上，通过matlab软件的功能，研究其数学物理方程的数值解，并通过对精确解和数值解进行对比，追究其数值解的可行性，在此，给出相关例子和程序代码，利于以后的再次研究和直接使用。[关键字]偏微分方程数值解matlab数学物理方程的可视化一：研究意义在我们解数学物理方程，理论上求数学物理方程的定解有着多种解法，但是有许多定解问题却不能严格求解，只能用数值方法求出满足实际需要的近似解。而且实际问题往往很复杂，这时即便要解出精确解就很困难，有时甚至不可能，另一方面，在建立数学模型时，我们已作了很多近似，所以求出的精确解也知识推导出的数学问题的精确解，并非真正实际问题的精确解。因此，我们有必要研究近似解法，只要使所求得的近似解与精确解之间的误差在规定的范围内，则仍能满足实际的需要，有限差分法和有限元法是两种最常用的求解数学物理方程的数值解法,而MATLAB在这一方面具有超强的数学功能，可以用来求其解。二：数值解法思想和步骤2.1:网格剖分为了用差分方法求解上述问题，将求解区域(,)|01,01xtxt作剖分。将空间区间[0,1]作m等分，将时间[0,1]区间作n等分，并记1/,1/,,0,,0jkhmnxjhjmtkkn。分别称h和为空间和时间步长。用两簇平行直线,0,,0jkxxjmttkn将分割成矩形网格。2.2：差分格式的建立0uutx………………………………（1）设G是,xt平面任一有界域，据Green公式(参考数学物理方程第五章)：()()Guudxdtudtudxtx其中G。于是可将（1）式写成积分守恒形式：()0udtudx…………………………（2）我们先从（2）式出发构造熟知的Lax格式设网格如下图所示图1取G为以(1,)Ajk，(1,1)Bjk，(1,1)Cjk和(1,)Djk为顶点的矩形。TABCDA为其边界，则()()()()()DABCABCDudtudxudxudxudtudt…………（3）右端第一个积分用梯形公式，第二个积分用中矩形公式，第三、四个积分用下矩形公式，则由（2）（3）式得到Lax-Friedrich格式111111()202kkkkkjjjjjuuuuuh截断误差为()[]kkkjhjjRuLuLu111111()22kkkkkkkjjjjjjjuuuuuuuhtx232223(),(0,0)26kkjjuuhOhjmkntx所以Lax-Friedrich格式的截断误差的阶式2()Oh令/sh：则可得差分格式为111111(),(0,0)222kkkkkjjjjjssuuuuujmkn0cos()(0)jjuxjmDCBAk+1kj+1j-1j0cos(),cos(),(0)kkkmkututkn2.3差分格式的求解差分格式111111(),(0,0)222kkkkkjjjjjssuuuuujmkn写成如下矩阵形式：10111212111100000222211000002222110000002222000000000000000kkkkkkmmkmkmssuussuuuussuu则需要通过对k时间层进行矩阵作用求出k+1时间层。对上面的矩阵形式我通过C++（或matlab）编出如附录的程序求出数值解、真实解和误差。例1：如下方程0,01,01.(,0)cos(),01.(0,)(1,)cos(),01.uuxttxuxxxututtt利用matlab的数值解法结果并填入表中关系对比如下：表11/900,1/1000(0.9hs）jk（x,t）数值解真实解误差450100(0.5,0.1)-0.308981-0.3090170.000036450200(0.5,0.2)-0.587647-0.5877850.000138450300(0.5,0.3)-0.808731-0.8090170.000286450400(0.5,0.4)-0.950609-0.9510560.000448450500(0.5,0.5)-0.999409-1.0000000.000591450600(0.5,0.6)-0.950496-0.9510570.000560450700(0.5,0.7)-0.808539-0.8090170.000478450800(0.5,0.8)-0.587437-0.5877850.000348450900(0.5,0.9)-0.308833-0.3090170.0001844501000(0.5,1.0)0.000002-0.0000000.000002从上面可以看出，数值解精度相当高的三：matlab的在数学物理方程上简单的应用Matlab是一个强大的计算工具，超强的计算能力和绘图能力，下面几个例题来说明matlab数学物理上的应用例1：将函数1/（1－a）2在z=0处展成幂级数。解：symsa；Taylor(1/(1-a)^2,0)ans=1+2*a+3*a^2+4*a^3+5*a^4+6*a^5例2：写出函数f（x）＝1/(x2＋p2)（a0）的Fourier变换式。解：symsxw;symsapotitivef=1/(x^2+p^2);F=Fourier(f,x,w)F=pi*(signum(0,Re(p),0)*cosh(p*w)-2*Heaviside(w)*sinh(p*w)+sinh(p*w))/p例2：已知函数f（x）＝x3e－x，试求取该函数的Laplace变换,并对结果进行Laplace反变换。解：symsxw;f=x^3*exp(-x);F=laplace(f,x,w)F=6/(w+1)^4对得出的结果进行Laplace反变换,从而有ilaplace(F)ans=x^3*exp(-x)利用手工方法对函数进行Fourier变换和Laplace变换，计算起来繁琐、复杂，且容易出错，利用MATLAB快速、准确。四：matlab解数学物理方程4.1：数值解法与精确解的可视化对比分析以下面的问题为例子（课本原题）根据上面的可建立方程如下：根据分离变量和差分其图形结果如下：图2分离变量热传导图3差分热传导其代码如下：图2x=0:0.1*pi:pi;y=0:0.4:10;[x,t]=meshgrid(x,y);u=0;m=length(j);%matlab可计算的最大数，相当于无穷fori=0:mu=u+8*(-1)^i/(pi*(2*i+1)^2)*(sin((2*i+1)/2*x).*exp(-(2*i+1)^2/4*t));end;surf(x,t,u);xlabel('x'),ylabel('t'),zlabel('T');title('分离变量法（无穷）');disp(u);图3u=zeros(20,100);%t=1x=pi20行100列横坐标为x纵坐标为ts=(1/100)/(pi/20)^2;fprintf('稳定性系数S为:\n');disp(s);fori=1:20u(i,1)=i/20*pi;;end;forj=1:100u(1,j)=0;endforj=1:99fori=2:19u(i,j+1)=s*u(i+1,j)+(1-2*s)*u(i,j)+s*u(i-1,j);endendforj=1:100u(20,j)=u(19,j);end;disp(u);[x,t]=meshgrid(1:100,1:20);surf(x,t,u);xlabel('t'),ylabel('x'),zlabel('T');title('有限差分法解');从上面可以看出数值解法精度很高，图形基本完全一样的4.2：matlab实现数值解法以下面的方程为例基本步骤：(1)区域的离散或子区域的划分；(2)插值函数的选择；(3)方程组的建立；(4)方程组的求解。a=input('请输入系数a的值：');l=input('请输入长度l的值：');M=input('请输入将区间[0,l]等分的个数M：');ot=input('请输入时间增量ot的值：');n=input('请输入运行次数n的值：');ox=l/M;x0=zeros(M+1,1);forii=1:Mx0(ii+1)=ii*ox;endu=sin(pi*x0/l);%t=0时u(x,t)的值r=a^2*ot/(ox)^2;forii=1:n%数据的输入B=zeros(M-1,1);%存放系数矩阵主对角线元素A=zeros(M-2,1);%存放系数矩阵主对角线元素下方次对角线的元素C=zeros(M-2,1);%存放系数矩阵主对角线元素上方次对角线的元素S=zeros(M-1,1);%存放右端的常数项forii=1:M-2B(ii)=1+2*r;A(ii)=-r;C(ii)=-r;S(ii)=u(ii+1,1);endB(M-1)=1+2*r;S(M-1)=u(M,1);u(1,2)=0;u(M+1,2)=0;S(1,1)=S(1,1)+r*u(1,2);S(M-1,1)=S(M-1,1)+r*u(M+1,2);%追赶法S(1)=S(1)/B(1);T=B(1);k=2;whilek~=MB(k-1)=C(k-1)/T;T=B(k)-A(k-1)*B(k-1);S(k)=(S(k)-A(k-1)*S(k-1))/T;k=k+1;endk=1;whilek~=M-1S(M-1-k)=S(M-1-k)-B(M-1-k)*S(M-k);k=k+1;endu(2:M,2)=S;%把结果放入矩阵u中u(:,1)=u(:,2);%过河拆桥end%计算精确值，存放在u的第二列forx=0:Mu(x+1,2)=exp(-(pi*a/l)^2*n*ot)*sin(pi*x*ox/l);end%计算最大相对误差ez=zeros(M-1,1);forii=2:Mez(ii-1)=abs(u(ii,1)-u(ii,2))/u(ii,2);endE=max(ez);fprintf('最后时刻数值解与精确解分别为:\n');disp(u);fprintf('差分法得到的结果与正确结果的最大相对误差为：');disp([num2str(E*100)'%']);%画二维图比较plot(x0,u(:,1),'r:',x0,u(:,2),'b-');legend('数值解','精确解')xlabel('x'),ylabel('u(x,t)')title('最后时刻热传导问题数值解与精确解比较')。运行如下：请输入系数a的值：1请输入长度l的值：1请输入将区间[0,l]等分的个数M：10请输入时间增量ot的值：0.001请输入运行次数n的值：100最后时刻数值解与精确解分别为:000.11670.11520.22190.21910.30540.30150.35910.35450.37750.37270.35910.35450.30540.30150.22190.21910.11670.115200.0000差分法得到的结果与正确结果的最大相对误差为：1.2934%图4数值解与精确解对比五:matlab偏微分工具箱介绍上面我们用的是编程技术，和c++没有多大区别，但是在matlab中还提供了偏微分工具箱，可以模拟数理方程动态图