数理方程数值解法与其在matlab软件上的实现张体强1026222 廖荣发1026226[摘要]数学物理方程的数值解在实际生活中越来越使用，首先基于偏微分数值解的思想上，通过matlab软件的功能，研究其数学物理方程的数值解，并通过对精确解和数值解进行对比，追究其数值解的可行性，在此，给出相关例子和程序代码，利于以后的再次研究和直接使用。

[关键字]偏微分方程数值解matlab 数学物理方程的可视化一：研究意义在我们解数学物理方程，理论上求数学物理方程的定解有着多种解法，但是有许多定解问题却不能严格求解，只能用数值方法求出满足实际需要的近似解。

而且实际问题往往很复杂，这时即便要解出精确解就很困难，有时甚至不可能，另一方面，在建立数学模型时，我们已作了很多近似，所以求出的精确解也知识推导出的数学问题的精确解，并非真正实际问题的精确解。

因此，我们有必要研究近似解法，只要使所求得的近似解与精确解之间的误差在规定的范围内，则仍能满足实际的需要，有限差分法和有限元法是两种最常用的求解数学物理方程的数值解法,而MATLAB 在这一方面具有超强的数学功能，可以用来求其解。

二：数值解法思想和步骤2.1:网格剖分为了用差分方法求解上述问题，将求解区域{}(,)|01,01x t x t Ω=≤≤≤≤作剖分。

将空间区间[0,1]作m 等分，将时间[0,1]区间作n 等分，并记1/,1/,,0,,0j k h m n x jh j m t k k n ττ===≤≤=≤≤。

分别称h 和τ为空间和时间步长。

用两簇平行直线,0,,0j k x x j m t t k n =≤≤=≤≤将Ω分割成矩形网格。

2.2：差分格式的建立0u u t x∂∂-=∂∂………………………………（1） 设G 是,x t 平面任一有界域，据Green 公式(参考数学物理方程第五章)：()()Gu udxdt udt udx t xΓ∂∂-=--∂∂⎰⎰ 其中G Γ=∂。

于是可将（1）式写成积分守恒形式：()0udt udx Γ--=⎰ (2)我们先从（2）式出发构造熟知的Lax 格式设网格如下图所示图1取G 为以(1,)A j k +，(1,1)B j k ++，(1,1)C j k -+和(1,)D j k -为顶点的矩形。

T ABCD = A 为其边界，则()()()()()DABCABCDudt udx u dx u dx u dt u dt Γ--=-+-+-+-⎰⎰⎰⎰⎰ (3)右端第一个积分用梯形公式，第二个积分用中矩形公式，第三、四个积分用下矩形公式，则由（2）（3）式得到Lax-Friedrich 格式111111()202k k k k k j j j j j u u u u u hτ+-+-+-+-+=截断误差为()[]k k kj h j j R u L u Lu =-111111()22k k k k k k k j j j j j j j u u u u u u u h t xτ+-+-+-+-∂∂=+-+∂∂232223(),(0,0)26kkjj u u h O h j m k n t xττ∂∂=-=+≤≤≤≤∂∂ 所以Lax-Friedrich 格式的截断误差的阶式2()O h τ+ 令/s h τ=：则可得差分格式为111111(),(0,0)222k k k kk j j j j j s s u u u u u j m k n +--++=-+++≤≤≤≤ 0cos()(0)j j u x j m π=≤≤k+10cos(),cos(),(0)k k k m k u t u t k n ππ==-≤≤2.3差分格式的求解差分格式111111(),(0,0)222k k k k k j j j j j s s u u u u u j m k n +--++=-+++≤≤≤≤写成如下矩阵形式：1011121211110000022221100022221100000022220000000000000k k k k k k m m k m k m s s u u s s u u u us s u u +++---⎛⎫-+ ⎪⎛⎫⎛⎫ ⎪ ⎪ ⎪ ⎪-+ ⎪ ⎪ ⎪⎪ ⎪ ⎪= ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪-+ ⎪ ⎪ ⎪ ⎪⎪ ⎪ ⎪⎝⎭⎝⎭ ⎪ ⎪⎝⎭则需要通过对k 时间层进行矩阵作用求出k+1时间层。

对上面的矩阵形式我通过C++（或matlab ）编出如附录的程序求出数值解、真实解和误差。

例1：如下方程0,01,0 1.(,0)cos(),0 1.(0,)(1,)cos(),0 1.u ux t t xu x x x u t u t t t ππ∂∂-=≤≤≤≤∂∂=≤≤=-=≤≤ 利用 matlab 的数值解法结果并填入表中关系对比如下：表11/900,1/1000(0.9h s τ===）从上面可以看出，数值解精度相当高的三：matlab的在数学物理方程上简单的应用Matlab是一个强大的计算工具，超强的计算能力和绘图能力，下面几个例题来说明matlab数学物理上的应用例1：将函数1/（1－a）2在z=0 处展成幂级数。

解：>>syms a；>>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 )（a>0）的Fourier 变换式。

解：>>syms x w;>>syms a potitive>>f=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）＝x3 e－x，试求取该函数的Laplace 变换,并对结果进行Laplace 反变换。

解：>>syms x w;>>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可计算的最大数，相当于无穷for i=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=1 x=pi 20行100列横坐标为x 纵坐标为t s=(1/100)/(pi/20)^2;fprintf('稳定性系数S为:\n');disp(s);for i=1:20u(i,1)=i/20*pi;;end;for j=1:100u(1,j)=0;endfor j=1:99for i=2:19u(i,j+1)=s*u(i+1,j)+(1-2*s)*u(i,j)+s*u(i-1,j);endendfor j=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);for ii=1:Mx0(ii+1)=ii*ox;endu=sin(pi*x0/l);%t=0 时u(x,t)的值r=a^2*ot/(ox)^2;for ii=1:n%数据的输入B=zeros(M-1,1);%存放系数矩阵主对角线元素A=zeros (M-2,1);%存放系数矩阵主对角线元素下方次对角线的元素C=zeros (M-2,1);%存放系数矩阵主对角线元素上方次对角线的元素S=zeros(M-1,1);%存放右端的常数项for ii=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;while k~=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;while k~=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 的第二列for x=0:Mu(x+1,2)=exp(-(pi*a/l)^2*n*ot)*sin(pi*x*ox/l); end%计算最大相对误差ez=zeros(M-1,1);for ii=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(' 最后时刻热传导问题数值解与精确解比较')。