数理方程基于matlab的数值解法.docx
《数理方程基于matlab的数值解法.docx》由会员分享，可在线阅读，更多相关《数理方程基于matlab的数值解法.docx（17页珍藏版）》请在冰豆网上搜索。
数理方程基于matlab的数值解法
数理方程数值解法与其在matlab软件上的实现
张体强1026222廖荣发1026226
[摘要]
数学物理方程的数值解在实际生活中越来越使用，首先基于偏微分数值解的思想上，通过matlab软件的功能，研究其数学物理方程的数值解，并通过对精确解和数值解进行对比，追究其数值解的可行性，在此，给出相关例子和程序代码，利于以后的再次研究和直接使用。
[关键字]
偏微分方程数值解matlab数学物理方程的可视化
一：
研究意义
在我们解数学物理方程，理论上求数学物理方程的定解有着多种解法，但是有许多定解问题却不能严格求解，只能用数值方法求出满足实际需要的近似解。
而且实际问题往往很复杂，这时即便要解出精确解就很困难，有时甚至不可能，另一方面，在建立数学模型时，我们已作了很多近似，所以求出的精确解也知识推导出的数学问题的精确解，并非真正实际问题的精确解。
因此，我们有必要研究近似解法，只要使所求得的近似解与精确解之间的误差在规定的范围内，则仍能满足实际的需要，有限差分法和有限元法是两种最常用的求解数学物理方程的数值解法,而MATLAB在这一方面具有超强的数学功能，可以用来求其解。
二：
数值解法思想和步骤
2.1:
网格剖分
为了用差分方法求解上述问题，将求解区域
作剖分。
将空间区间
作
等分，将时间
区间作
等分，并记
。
分别称
和
为空间和时间步长。
用两簇平行直线
将
分割成矩形网格。
2.2：
差分格式的建立
………………………………
（1）
设
是
平面任一有界域，据Green公式（参考数学物理方程第五章）：
其中
。
于是可将
（1）式写成积分守恒形式：
…………………………
（2）
我们先从
（2）式出发构造熟知的Lax格式设网格如下图所示
图1
取
为以
，
，
和
为顶点的矩形。
为其边界，则
…………（3）
右端第一个积分用梯形公式，第二个积分用中矩形公式，第三、四个积分用下矩形公式，则由
（2）（3）式得到Lax-Friedrich格式
截断误差为
所以Lax-Friedrich格式的截断误差的阶式
令
：
则可得差分格式为
2.3差分格式的求解
差分格式
写成如下矩阵形式：
则需要通过对k时间层进行矩阵作用求出k+1时间层。
对上面的矩阵形式我通过C++（或matlab）编出如附录的程序求出数值解、真实解和误差。
例1：
如下方程
利用matlab的数值解法结果并填入表中关系对比如下：
表1
j
k
（x,t）
数值解
真实解
误差
450
100
（0.5,0.1）
-0.308981
-0.309017
0.000036
450
200
（0.5,0.2）
-0.587647
-0.587785
0.000138
450
300
（0.5,0.3）
-0.808731
-0.809017
0.000286
450
400
（0.5,0.4）
-0.950609
-0.951056
0.000448
450
500
（0.5,0.5）
-0.999409
-1.000000
0.000591
450
600
（0.5,0.6）
-0.950496
-0.951057
0.000560
450
700
（0.5,0.7）
-0.808539
-0.809017
0.000478
450
800
（0.5,0.8）
-0.587437
-0.587785
0.000348
450
900
（0.5,0.9）
-0.308833
-0.309017
0.000184
450
1000
（0.5,1.0）
0.000002
-0.000000
0.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）（a>0）的Fourier变换式。
解：
>>symsxw;
>>symsapotitive
>>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）＝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差分热传导
其代码如下：
图2
x=0:
0.1*pi:
pi;
y=0:
0.4:
10;
[x,t]=meshgrid（x,y）;
u=0;
m=length（j）;%matlab可计算的最大数，相当于无穷
fori=0:
m
u=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）;
图3
u=zeros（20,100）;%t=1x=pi20行100列横坐标为x纵坐标为t
s=（1/100）/（pi/20）^2;
fprintf（'稳定性系数S为:
\n'）;
disp（s）;
fori=1:
20
u（i,1）=i/20*pi;;
end;
forj=1:
100
u（1,j）=0;
end
forj=1:
99
fori=2:
19
u（i,j+1）=s*u（i+1,j）+（1-2*s）*u（i,j）+s*u（i-1,j）;
end
end
forj=1:
100
u（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:
M
x0（ii+1）=ii*ox;
end
u=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-2
B（ii）=1+2*r;A（ii）=-r;C（ii）=-r;
S（ii）=u（ii+1,1）;
end
B（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~=M
B（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;
end
k=1;
whilek~=M-1
S（M-1-k）=S（M-1-k）-B（M-1-k）*S（M-k）;
k=k+1;
end
u（2:
M,2）=S;%把结果放入矩阵u中
u（:
1）=u（:
2）;%过河拆桥
end
%计算精确值，存放在u的第二列
forx=0:
M
u（x+1,2）=exp（-（pi*a/l）^2*n*ot）*sin（pi*x*ox/l）;
end
%计算最大相对误差
ez=zeros（M-1,1）;
forii=2:
M
ez（ii-1）=abs（u（ii,1）-u（ii,2））/u（ii,2）;
end
E=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
最后时刻数值解与精确解分别为:
00
0.11670.1152
0.22190.2191
0.30540.3015
0.35910.3545
0.37750.3727
0.35910.3545
0.30540.3015
0.22190.2191
0.11670.1152
00.0000
差分法得到的结果与正确结果的最大相对误差为：
1.2934%
图4数值解与精确解对比
五:
matlab偏微分工具箱介绍
上面我们用的是编程技术，和c++没有多大区别，但是在matlab中还提供了偏微分工具箱，可以模拟数理方程动态图，这里简单介绍一些功能
5.1：
偏微分方程工具箱的功能
偏微分方程工具箱（PDEToolbox）提供了研究和求解空间二维偏微分方程问题的一个强大而又灵活实用的环境。
PDEToolbox的功能包括：
（1）设置PDE（偏微分方程）定解问题，即设置二维定解区域、边界条件以及方程的形式和系数；
（2）用有限元法（FEM）求解PDE数值解；
（3）解的可视化。
下面我简单介绍一下操作
启动matlab，在工作区输入pdetool，会出现一个界面，如下
图5PDEToolbox界面
系统立即产生偏微分方程工具箱（PDEToolbox）的图形用户界面（GraphicalUserInterface，简记为GUI），即PDE解的图形环境，这时就可以在它上面画出定解区域、设置方程和边界条件、作网格剖分、求解、作图等工作。
5.2：
偏微分工具箱的求解范围
5.2.1：
方程类型
PDEToolbox求解的基本方程有椭圆型方程、抛物型方程、双曲型方程、特征值方程、椭圆型方程组以及非线性椭圆型方程。
（参考这学期数学物理方程几大类型）另外，对于Poission方程还有一个矩形网格的快速求解器
5.2.2如何使用FDEToolbox
1：
定解问题的设置
员简单的办法是在PDETool上直接使用图形用户界面（GUl）。
设置定解问题包括三个步骤：
（1）Draw模式：
使用CSG（几何结构实体模型）对话框画几何区域，包括矩形、圆、椭圆和多边形，也可以将它们组合使用。
（2）Boundary模式：
在各个边界段上给出边界条件，
（3）PDE模式：
确定方程的类型、系数c，a，f和dc。
也能够在不同子区域上设置不同的系数（反映材料的性质）。
5.2.3解PDE问题
用GUI解PDE问题主要经过下面两个过程（模式）
（1）Mesh模式；生成网格．自动控制网格参数。
（2）Solve模式：
对于椭圆型方程还能求非线性和自适应解。
对于抛物型和双曲型力程．设置初始边值条件后能求出给定t时刻的解。
对于特征值问题，能求出给定区间内的特征值；求解后可以加密网格再求解。
5.2.3使用Toolbox求解非标准的问题
对于非标准的问题。
可以用PDEToo1box的函数。
或者用FEM（有限元法）求解更为复杂的问题。
5.2.4计算结果的可视化
从GUI能够使用Plot模式实现可视化。
可以使用Color,Height和Vector等作图。
对于抛物型和双曲型方程，还可以生成解的动画。
这些操作通过命令行都很容易实现。
5.2.5应用领域
在应用界面提供了丁如下应用领域
．结构力学——平面应力问题
．结构力学——平面应变问题
．静电场问题
．静磁场问题
．交流电磁场问题
．直流导体介质问题
．热传导问题
下面是一个我们模拟一个简单的二维波动方程，其中程序部分如下
[p,e,t]=initmesh（‘squareg’）;
x=p（1,:
）’;y=p（2,:
）’;
u0=atan（cos（pi/2*x））;
ut0=3*sin（pi*x）.*exp（sin（pi/2*y））;
n=3l;
tlist=linspace（0,5,n）;
uu=hyperbolic（u0,ut0,tlist,‘squareb3’,p,e,t,1,0,0,1）;
delta=-1:
0.1:
1;
[uxy,tn,a2,a3]=tri2grid（p,t,uu（:
1）,delta,delta）;
gp=[tn;a2;a3];
umax=max（max（uu））;
umin=min（min（uu））;
newplot
M=moviein（n）;
ForI=1:
n,
pdeplot（p,e,t,‘xydata’,uu（:
I）,‘zdata’,uu（:
I）,‘mesh’,‘off’,‘xygrid’,‘on’,‘gridparam’,gp,‘colorbar’,‘off’,‘zstyle’,‘continuous’）;
axis（[-11–11uminumax]）;caxis（[uminumax]）;
M（:
I）=getframe;
End
Movie（M,10）;
图6matlab模拟二维波动方程
六：
总结
在以后我们现实生活中很少出现非数值解，数值解法的思想很重要，在基于思想上我们可以通过计算机实现，matlab在这一方面非常方便，可以通过纯的编程实现其解的值，也可以借助其工具箱解方程，可以模拟二维图，三维图和动态图，学习不但只是学，还要懂得和生活练习，去实践，而本文章通过课本上的例题进而分析数值解法与精确解的区别，从而认为数值解不但可行，而且很使用，在计算机实现数学物理方程的数值解时，只要对以上的程序稍作修改，即可得数值解，还可以通过工具箱模拟出来，为科学研究提供一个更好的平台。
参考文献
[1]曹钢,王桂珍,任晓荣.一维热传导方程的基本解[J].山东轻
工业学院学报,2005,19（4）:
76-80.
[2]万正苏，方春华，张再云.关于热传导方程有限差分区域分
解并行算法精度的注记[J].湖南理工学院学报（自然科学
版），2007,20（3）：
12-14.
[3]StephenJ.Chapman.MATLAB编程[M].邢树军，郑碧波，译.
北京：
科学出版社，2008.
[4]田兵.用MATLAB解偏微分方程[J].阴山学刊，2006,20（4）:
12-13.
[5]王飞，裴永祥.有限差分方法的MATLAB编程[J].新疆师
范大学学报（自然科学版），2003,22（4）:
21-27.
[6]王宝红.热传导方程的可视化探讨[J].忻州师范学院学报,
2008,24
（2）:
31-36.
[7]李先枝.热传导方程差分解法的最佳网格[J].河南大学学报
（自然科学版），2004,34（3）:
16-18.
[8]赵德奎，刘勇.MATLAB在有限差分数值计算中的应用[J].
四川理工学院学报，2005，18（4）:
61-64.
[9]谢焕田，吴艳.拉普拉斯有限差分法的MATLAB实现[J].
四川理工学院学报，2008，21（3）:
1-2.
[10]赵刚，严尚安等主编，《数学物理方程》，中国科技出版社，2004年版，第8页
[11]复旦大学数学系主编，《数学物理方程》，上海科学技术出版社，1960年版，第10页
[12]闻新等编著，《MTLAB科学图形构建基础与应用》，科学出版社，2002年版，第107页
[13]ZakarauskasRierre,DossoStanleyE,FawcettJohnA.Matched-fieldinversionforsourcelocationandoptimaliquivalentbathymetry[J].J
[14]梁昆淼编，《数学物理方法》，人民教育出版社1978年第二版，第256—257页，第390—392页
[15]彭芳麟著，《数学物理方程的MATLAB解法与可视化》，清华大学出版社，2004年版，第140页，第142页，第150页
[16]李好,杨天春,王齐仁.基于Matlab7.0PDE工具箱求解数学物理方程[J].电脑开发与应用,2009,22
（1）:
26-27,48.
2012年5月
升级会员 