数理方程解的仿真.docx
《数理方程解的仿真.docx》由会员分享，可在线阅读，更多相关《数理方程解的仿真.docx（17页珍藏版）》请在冰豆网上搜索。
数理方程解的仿真
MATLAB实现
数理方程解的仿真
摘要
数学物理方程的数值解在实际生活中越来越使用，首先基于偏微分数值解的思想上，通过MATLAB软件的功能，研究其数学物理方程的数值解，并通过对精确解和数值解进行对比，追究其数值解的可行性。
MATLAB是高性能的数值计算型数学类科技应用软件，具有优秀的数值计算功能和强大的数据可视化能力。
在《数学物理方法》中应用MATLAB进行习题求解和计算机仿真，一方面可以提高解题的速度，另一方面可将抽象的解和一些特殊函数以图形形式显示出来，直观明了，物理意义明确。
关键词：
数学物理方法，可视化，MATLAB，图形绘制
引言
物理问题的解决最终往往是一个或若干个数学问题，如代数方程（组）、函数方程、微分方程、定积分、特殊函数、矩阵本征值、线性拟合、插值、傅里叶变换等。
传统物理学包括实验物理和理论物理两个部分。
实验物理以实验和观测为手段来发现新的物理现象，为理论物理提供新的理论规律的素材，并检验理论物理推论的正确程度及应用范围等等。
而理论物理是则从一系列的基本物理原理出发，列出数学方程，并用传统的数学分析方法求出显式的解析解，利用这些显式解析解所得出的结论与实验和观测结果对比，然而很多实际问题用传统的物理学既无法得到其解析解，又不能通过实验的办法取得经验和数据。
这些问题的出现不得不让我们思考该怎么办。
在我们解数学物理方程，理论上求数学物理方程的定解有着多种解法，但是有许多定解问题却不能严格求解，只能用数值方法求出满足实际需要的近似解。
而且实际问题往往很复杂，这时即便要解出精确解就很困难，有时甚至不可能，另一方面，在建立数学模型时，我们已作了很多近似，所以求出的精确解也知识推导出的数学问题的精确解，并非真正实际问题的精确解。
因此，我们有必要研究近似解法，只要使所求得的近似解与精确解之间的误差在规定的范围内，则仍能满足实际的需要，有限差分法和有限元法是两种最常用的求解数学物理方程的数值解法,而MATLAB在这一方面具有超强的数学功能，可以用来求其解。
矩阵是指纵横排列的二维数据表格，最早来自于方程组的系数及常熟构成的方阵。
MATLAB——矩阵实验Matrix Laboratory的简称，MATLAB的基本数据单位是矩阵，主要用于算法开发看，数据可视化，数据分析以及数值计算的高级技术计算语言和交互式环境，主要包括MATLAB和Simulink两大部分。
它将数值分析、矩阵计算、科学数据可视化以及非线性动态系统的建模和仿真等诸多强大功能集成在一个易于使用的视窗环境中，为科学研究、工程设计以及必须进行有效数值计算的众多科学领域提供了一种全面的解决方案。
第一章数值解原理
第1节网格剖分
为了用差分方法求解上述问题，将求解区域
作剖分。
将空间区间
作
等分，将时间
区间作
等分，并记
。
分别称
和
为空间和时间步长。
用两簇平行直线
将
分割成矩形网格。
第2节差分格式的建立
………………………………
（1）
设
是
平面任一有界域，据Green公式（参考数学物理方程第五章）：
其中
。
于是可将
（1）式写成积分守恒形式：
…………………………
（2）
我们先从
（2）式出发构造熟知的Lax格式设网格如下图所示
图1
取
为以
，
，
和
为顶点的矩形。
T=ABCD,
为其边界，则
…………（3）
右端第一个积分用梯形公式，第二个积分用中矩形公式，第三、四个积分用下矩形公式，则由
（2）（3）式得到Lax-Friedrich格式
截断误差为
所以Lax-Friedrich格式的截断误差的阶式
令
：
则可得差分格式为
第3节差分格式的求解
差分格式
写成如下矩阵形式：
则需要通过对k时间层进行矩阵作用求出k+1时间层。
第二章MATLAB仿真与可视化
第1节有界弦振动问题解析解（精确解）的图形分布
（1）
解析解：
（取系数为1）
（2）
1）n=1本振解随时间变化的图形分布
MATLAB仿真程序：
clear
a=1;l=1;x=0:
0.05:
1;t=0:
0.001:
3;u=0;n=1;%函数参数
[X,T]=meshgrid（x,t）;%函数变量
u=（cos（n*pi*a*T/l）+sin（n*pi*a*T/l））.*sin（n*pi*X/l）%函数定义式
figure
（1）
axis（[01-0.050.05]）%图形取点，x分布
mesh（X,T,u）
title（'n=1本振随时分布的透视图'）%图像名
xlabel（'x'）%图像坐标轴'x'
ylabel（'t'）%图像坐标轴't'
zlabel（'u'）%图像坐标轴'u'
MATLAB可视化仿真结果：
图2n=1本振解随时间变化的图形
2）n=2本振解随时间变化的图形分布
MATLAB仿真程序：
clear
a=1;l=1;x=0:
0.05:
1;t=0:
0.001:
3;u=0;n=2;%函数参数
[X,T]=meshgrid（x,t）;%函数变量¿
u=（cos（n*pi*a*T/l）+sin（n*pi*a*T/l））.*sin（n*pi*X/l）%函数定义式
figure
（1）
axis（[01-0.050.05]）%图形取点，x分布
mesh（X,T,u）
title（'n=2本振随时分布的透视图'）%图像名
xlabel（'x'）%图像坐标轴'x'
ylabel（'t'）%图像坐标轴'y'
zlabel（'u'）%图像坐标轴'u'
MATLAB可视化仿真结果：
图3n=2本振解随时间变化的图形
第2节数值解法与精确解的可视化对比分析
以下面的问题为例子
根据上面的可建立方程如下：
解析解的MATLAB源代码：
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）;
MATLAB可视化仿真结果：
图4分离变量热传导
数值解的MATLAB源代码：
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（'有限差分法解'）;
MATLAB可视化仿真结果：
图5差分热传导
数值解是采用某种计算方法,如有限元的方法,数值逼近,插值的方法,得到的解。
解析解就是一些严格的公式，给出任意的自变量就可以求出因变量，也称作精确解。
而上面差分法则采用了有限变量的数值解，分离变量法则采用无限变量的精确解，观察可视化图形结果可以得出数值解法精度很高，图形基本完全一致。
第3节MATLAB实现数值解法
以下面的方程为例
基本步骤：
（1）区域的离散或子区域的划分；
（2）插值函数的选择；
（3）方程组的建立；
（4）方程组的求解。
MATLAB源代码：
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（'最后时刻热传导问题数值解与精确解比较'）
MATLAB仿真可视化结果：
图6数值解与精确解对比
程序解的运行如下：
>>shuzhijiefa
请输入系数a的值：
1
请输入长度l的值：
1
请输入将区间[0,l]等分的个数M：
10
请输入时间增量ot的值：
0.0014
请输入运行次数n的值：
100
最后时刻数值解与精确解分别为:
00
0.07920.0776
0.15070.1476
0.20740.2032
0.24380.2388
0.25640.2511
0.24380.2388
0.20740.2032
0.15070.1476
0.07920.0776
00.0000
差分法得到的结果与正确结果的最大相对误差为：
2.0847%
第三章总结与体会
通过以上各项可视化以及代数借的仿真比对，可以明确得出，数值解与精确解很相似，而利用MATLAB工具有可以很好地实现数值解的仿真，所以将MATLAB应用在数理方程的求解上，是一种很重要的方法。
数值解法的思想很重要，在基于思想上我们可以通过计算机实现，MATLAB在这一方面非常方便，可以通过纯的编程实现其解的值，也可以借助其工具箱解方程，可以模拟二维图，三维图和动态图，学习不但只是学，还要懂得和生活练习，去实践，而本文章通过课本上的例题进而分析数值解法与精确解的区别，从而认为数值解不但可行，而且很实用，在计算机实现数学物理方程的数值解时，只要对以上的程序稍作修改，即可得数值解，还可以通过工具箱模拟出来，为科学研究提供一个更好的平台。
升级会员 