微分方程数值解实验.docx
《微分方程数值解实验.docx》由会员分享，可在线阅读，更多相关《微分方程数值解实验.docx（22页珍藏版）》请在冰点文库上搜索。
微分方程数值解实验
微分方程数值解
课程设计报告
班级：
______________
姓名：
_________
学号：
___________
成绩：
2017年6月21日
一、摘要1
二、常微分方程数值解2
2.14阶Runge-Kutta法和Adams4阶外插法的基本思路2
2.2算法流程图2
2.3用matlab编写源程序2
2.4常微分方程数值解法应用举例4
三、常系数扩散方程的经典差分格式6
3.1有限差分法的基本思路6
3.2算法流程图7
3.3用matlab编写源程序7
3.4有限差分法应用举例8
四、椭圆型方程的五点差分格式10
4.1五点差分法的基本思路10
4.2算法流程图11
4.3用matlab编写源程序11
4.4五点差分法应用举例12
五、自我总结16
六、参考文献16
一、摘要
自然界与工程技术中的很多现象，可以归结为微分方程定解问题。
其中，常微分方程求解是微分方程的重要基础内容。
但是，对于许多的微分方程，往往很难得到甚至不存在精确的解析表达式，这时候，数值解提供了一个很好的解决思路。
，针对于此，本文对常微分方程数值解法进行了简单研究，主要讨论了一些常用的数值解法，如欧拉法、改进的欧拉法、Runge—Kutta方法、Adams法以及椭圆型方程、抛物型方程的有限差分方法等，通过具体的算例，结合MATLAB求解画图，初步给出了一般常微分方程数值解法的求解过程。
同时，通过对各种方法的误差分析，让大家对各种方法的特点和适用范围有一个直观的感受。
关键词：
微分方程数值解、MATLAB
二、常微分方程数值解
2.1基本思路
常微分方程数值解法（numericalmethodsforordinarydifferentialequations）计算数学的一个分支.是解常微分方程各类定解问题的数值方法.现有的解析方法只能用于求解一些特殊类型的定解问题，实用上许多很有价值的常微分方程的解不能用初等函数来表示，常常需要求其数值解.所谓数值解，是指在求解区间内一系列离散点处给出真解的近似值.这就促成了数值方法的产生与发展.
常微分方程初值问题的数值解法是求方程
（1）的解在点列
上的近似值
，这里
是
到
的步长，一般略去下标记为
。
（1）
经典的
方法是一个四阶的方法，它的计算公式是：
（2）
方法的优点是：
单步法、精度高，计算过程便于改变步长，缺点是计算量较大，每前进一步需要计算四次函数值
。
在用龙格库塔方法时,要注意
的选择要合适，
太大，会使计算量加大，
太小，
较大，可能会使误差增大。
因此选择合适的
很重要。
我们要在考虑精度的基础上，选择合适的
。
2.2算法步骤
2.2.1、四阶龙格－库塔（R-K）方法流程图：
2.2.2、Adams4阶外插法流程图：
实例求解流程：
2.3用matlab编写源程序
Matlab程序源代码：
--------------定义Rk4.m文件----------------------------------------
functiondy=Rk4（x,y）
dy=zeros（3,1）;
dy
（1）=10*（-y
（1）+y
（2））;
dy
（2）=28*y
（1）-y
（2）-y
（1）*y（3）;
dy（3）=y
（1）*y
（2）-8*y（3）/3;
end
-------------------------------------------------------------------------------------------
--------------定义Adams4.m文件---------------------------------------
function[x,y,z]=adams4（x1,y1,z1,x2,y2,z2,x3,y3,z3,h）
%Adams外插法
kfy=0;ksy=0;kty=0;
kfz=0;ksz=0;ktz=0;
kfx=10*（y3-x3）;%eval_r（abx）;
kfy=x3*（28-z3）-y3;%eval_r（aby）;
kfz=x3*y3-8/3*z3;%eval_r（abz）;
ksx=10*（y2-x2）;
ksy=x2*（28-z2）-y2;;
ksz=x2*y2-8/3*z2;
ktx=10*（y1-x1）;
kty=x1*（28-z1）-y1;
ktz=x1*y1-8/3*z1;
x=x3+h/12*（23*kfx-16*ksx+5*ktx）;
y=y3+h/12*（23*kfy-16*ksy+5*kty）;
z=z3+h/12*（23*kfz-16*ksz+5*ktz）;
end
-------------------------------------------------------------------------------------------
--------------定义exe11.m文件------------------------------------------
[t,y]=ode45（Rk4,[0,30],[12,2,9]）
suptitle（'Runge-Kutta4阶法'）%总标题
subplot（2,2,1）;
plot（t,y（:
1））;gridon;
legend（'x关于t的变化关系图',1）;
xlabel（'t','FontSize',14）;
ylabel（'x','FontSize',14）;
subplot（2,2,2）;
plot（t,y（:
2））;gridon;
legend（'y关于t的变化关系图',1）;
xlabel（'t','FontSize',14）;
ylabel（'y','FontSize',14）;
subplot（2,2,3）
plot（t,y（:
3））;gridon;
legend（'z关于t的变化关系图',1）;
xlabel（'t','FontSize',14）;
ylabel（'z','FontSize',14）;
subplot（2,2,4）
plot3（y（:
1）,y（:
2）,y（:
3））;gridon;
legend（'x,y,z的空间关系图',1）;
xlabel（'x','FontSize',14）;
ylabel（'y','FontSize',14）;
zlabel（'y','FontSize',14）;
view（40,60）;%锁定同样的视图，便于比较
-------------------------------------------------------------------------------------------
--------------定义exe12.m文件------------------------------------------
a=1:
1:
3000;b=1:
1:
3000;c=1:
1:
3000;
t=1:
1:
3000;
x1=5;y1=5;z1=10;
x2=7.83;y2=14.30;z2=12.34;
x3=15.32;y3=22.87;z3=28.07;
fori=1:
1:
3000
[x,y,z]=adams4（x1,y1,z1,x2,y2,z2,x3,y3,z3,0.01）;
a（i）=x;b（i）=y;c（i）=z;
x1=x2;y1=y2;z1=z2;
x2=x3;y2=y3;z2=z3;
x3=x;y3=y;z3=z;
fprintf（'x=%f''y=%f''z=%f\n',x,y,z）;%显示迭代值
end
suptitle（'Adams4阶外插法'）%总标题
subplot（2,2,1）;
plot（t,a）;gridon;
legend（'x关于t的变化关系图',1）;
xlabel（'t','FontSize',14）;
ylabel（'x','FontSize',14）;
subplot（2,2,2）;
plot（t,b）;gridon;
legend（'y关于t的变化关系图',1）;
xlabel（'t','FontSize',14）;
ylabel（'y','FontSize',14）;
subplot（2,2,3）;
plot（t,c）;gridon;
legend（'z关于t的变化关系图',1）;
xlabel（'t','FontSize',14）;
ylabel（'z','FontSize',14）;
subplot（2,2,4）
plot3（a,b,c）;gridon;
legend（'x,y,z的空间关系图',1）;
xlabel（'x','FontSize',14）;
ylabel（'y','FontSize',14）;
zlabel（'y','FontSize',14）;
view（40,60）;%锁定同样的视图，便于比较
-------------------------------------------------------------------------------------------
2.4常微分方程数值解法应用举例
题目：
MIT的气象学家洛仑兹（E.Lorenz）在1963年研究大气对天气的影响时，提出了Lorenz方程：
其中
，
，
此方程初值问题的解存在且唯一。
当
时，Lorenz方城有两个不稳定的不动点，
和
，一个不稳定的不动点
。
当
时，
和
都变成不稳定的。
此时，存在混沌和一个奇怪吸引子。
解：
由于
时，
和
都变成不稳定的。
此时，存在混沌和一个奇怪吸引子。
所以我们取
来专门研究该现象。
利用Matlab程序计算：
命令窗口输入：
Runge-Kutta4阶法：
>>exe11
Adams4阶外插法：
>>exe12
结果输出：
由于迭代次数的关系，结果数据量很大，故这里以图片来展示结果。
Runge-Kutta4阶法结果：
Adams4阶外插法结果：
三、常系数扩散方程的经典差分格式
3.1有限差分法的基本思路
用有限差分法解常系数扩散方程
有加权隐式差分格式
其中
，当
时为Crank-Nicolson格式，当
时为向后差分格式，当
时为向前差分格式。
加权隐式格式稳定的条件是
加权隐式格式是两层隐式格式，用第n层计算第n+1层节点值的时候，要解线性方程组。
3.2算法步骤
实验步骤如下：
（1）输入
，确定加权隐式格式的参数；
（2）定义向量v，把初边值条件离散，得到
的值存入向量v；
（3）利用差分格式由第n层计算第n+1层，建立相应线性方程组，求解并且存入向量v；
（4）计算到
，输出
。
3.3用matlab编写源程序
Matlab程序源代码：
--------------定义parabola.m文件------------------------------------------
function[u,tt,uu]=parabola（k,q,l）%q=thetal=tao
h=0.1;%h:
空间步长
t=l*h*h;
x=0.1:
0.1:
0.9;%x坐标取值
forn=1:
9
u（n）=sin（x（n）*pi）;%初始值
end
u=u';
v=zeros（9,1）;%初始化v矩阵
forn=1:
k
A=zeros（9,9）;
A（1,1）=1+2*q*t/h/h;
fori=2:
9
A（i,i）=1+2*q*t/h/h;
A（i-1,i）=-q*t/h/h;
A（i,i-1）=-q*t/h/h;
end
B=zeros（9,1）;
B（1,1）=（1-q）*t/h/h*（u
（2）-2*u
（1））+u
（1）;
B（9,1）=（1-q）*t/h/h*（u（8）-2*u（9））+u（9）;
fori=2:
8
B（i,1）=（1-q）*t/h/h*（u（i-1）-2*u（i）+u（i+1））+u（i）;%加权隐式格式
end
v=inv（A）*B;%求解矩阵v
w=u;
u=v;
v=w;
n=n+1;
end
tt=k*t;
uu=zeros（9,1）;%初始化uu矩阵
fori=1:
9
uu（i,1）=sin（pi*x（i））*exp（-pi*pi*tt）;%精确值
end
-------------------------------------------------------------------------------------------
--------------定义exe20.m文件------------------------------------------
[u,tt,uu]=parabola（100,0.25,0.1）;
fprintf（'τ=0.25，θ=0.1时0.4处的准确值为：
%f\n',uu（4））
[u,tt,uu]=parabola（100,0.25,0.1）;
fprintf（'τ=0.25，θ=0.1时0.4处的数值值为：
%f\n',u（4））
fprintf（'误差为：
%f\n\n\n',abs（uu（4）-u（4）））
x=0.1:
0.1:
0.9;
subplot（2,2,1）;
plot（x,u,x,uu,'o'）;
gridon;
title（'τ=0.25，θ=0.1'）;
[u,tt,uu]=parabola（100,0.25,0.5）;
fprintf（'τ=0.25，θ=0.5时0.4处的准确值为：
%f\n',uu（4））
[u,tt,uu]=parabola（100,0.25,0.5）;
fprintf（'τ=0.25，θ=0.5时0.4处的数值值为：
%f\n',u（4））
fprintf（'误差为：
%f\n\n\n',abs（uu（4）-u（4）））
x=0.1:
0.1:
0.9;
subplot（2,2,2）;
plot（x,u,x,uu,'o'）;
gridon;
title（'τ=0.25，θ=0.5'）;
[u,tt,uu]=parabola（100,0.5,0.1）;
fprintf（'τ=0.5，θ=0.1时0.4处的准确值为：
%f\n',uu（4））
[u,tt,uu]=parabola（100,0.5,0.1）;
fprintf（'τ=0.5，θ=0.1时0.4处的数值值为：
%f\n',u（4））
fprintf（'误差为：
%f\n\n\n',abs（uu（4）-u（4）））
x=0.1:
0.1:
0.9;
subplot（2,2,3）;
plot（x,u,x,uu,'o'）;
gridon;
title（'τ=0.5，θ=0.1'）;
[u,tt,uu]=parabola（100,0.5,0.5）;
fprintf（'τ=0.5，θ=0.5时0.4处的准确值为：
%f\n',uu（4））
[u,tt,uu]=parabola（100,0.5,0.5）;
fprintf（'τ=0.5，θ=0.5时0.4处的数值值为：
%f\n',u（4））
fprintf（'误差为：
%f\n\n\n',abs（uu（4）-u（4）））
x=0.1:
0.1:
0.9;
subplot（2,2,4）;
plot（x,u,x,uu,'o'）;
gridon;
title（'τ=0.5，θ=0.5'）;
-------------------------------------------------------------------------------------------
3.4有限差分法应用举例
题目：
考虑常系数扩散方程的初边值问题
其解析解为
用加权隐式格式近似求解
其中
，取
为时间步长，
为网格比，对不同的时间步长（
），计算当
时初边值问题的解
，并且与精确解比较，分析比较结果。
用有限差分法近似求解，并且与精确解比较，分析结果。
解：
利用Matlab程序计算：
命令窗口输入：
exe20
结果输出：
τ=0.25，θ=0.1时0.4处的准确值为：
0.354466
τ=0.25，θ=0.1时0.4处的数值值为：
0.356486
误差为：
0.002020
τ=0.25，θ=0.5时0.4处的准确值为：
0.006840
τ=0.25，θ=0.5时0.4处的数值值为：
0.006696
误差为：
0.000143
τ=0.5，θ=0.1时0.4处的准确值为：
0.354466
τ=0.5，θ=0.1时0.4处的数值值为：
0.357343
误差为：
0.002877
τ=0.5，θ=0.5时0.4处的准确值为：
0.006840
τ=0.5，θ=0.5时0.4处的数值值为：
0.007115
误差为：
0.000275
误差
结论：
根据图和数据可以看出：
当
，
时，得到的结果误差最小。
当改变q值和k值，可得到相应的结果。
四、椭圆型方程的五点差分格式
4.1五点差分法的基本思路
对Laplace方程的第一边值问题
利用taylor展开可得逼近它的五点差分格式的差分逼近
其中
分别为
轴和
轴步长，边界条件可以由
离散可得，当
时有
。
注意五点格式计算节点是由边界的已知节点，计算内部节点，计算时需要联立大型方程组，该方程组可以用迭代法求解。
4.2算法步骤
主要步骤：
（1）首先取定
，对求解区域划分网格，按照网格定义矩阵v1，使得矩阵里面每个元素对应求解区域中的每个节点；
（2）由边界条件定义v1矩阵中边界元素的值，其余元素定义为零；
（3）定义与v1同型的零矩阵v2；
（4）用五点格式公式通过矩阵v1迭代计算矩阵v2，迭代精度为0.1；
（5）画图。
4.3用matlab编写源程序
Matlab程序源代码：
--------------定义fivepoint.m文件------------------------------------------
function[peuxyk]=fivepoint（h,m,n,kmax,ep）
%h步长
%m,n为x,y方向的网格数，例如（2-0）/0.01=200;
%kmax为最大迭代次数
%e为误差，p为精确解
symstemp;
u=zeros（n+1,m+1）;
x=0+（0:
m）*h;
y=0+（0:
n）*h;
for（i=1:
n+1）
u（i,1）=sin（pi*y（i））;%（0,y）
u（i,m+1）=exp
（1）*exp
（1）*sin（pi*y（i））;%（17,y）
end
for（i=1:
n）
for（j=1:
m）
f（i,j）=（pi*pi-1）*exp（x（j））*sin（pi*y（i））;%求区域中方程的每个节点
end
end
t=zeros（n-1,m-1）;
for（k=1:
kmax）
for（i=2:
n）
for（j=2:
m）%五点差分
temp=h*h*f（i,j）/4+（u（i,j+1）+u（i,j-1）+u（i+1,j）+u（i-1,j））/4;
t（i,j）=（temp-u（i,j））*（temp-u（i,j））;
u（i,j）=temp;
end
end
t（i,j）=sqrt（t（i,j））;
if（k>kmax）
break;%达到最大迭代次数，终止循环
end
if（max（max（t））break;
end
end
for（i=1:
n+1）
for（j=1:
m+1）
p（i,j）=exp（x（j））*sin（pi*y（i））;%精确解
e（i,j）=abs（u（i,j）-exp（x（j））*sin（pi*y（i）））;%误差
end
end
--------------------------------------------------------------------------------------------
4.4五点差分法应用举例
题目：
给定如下Laplace方程（poisson方程的特殊情况）的定解问题
利用椭圆型方程的五点格式，计算该问题的近似解，并且画出近似解的图形。
解：
利用Matlab程序计算：
命令窗口输入：
[peuxyk]=fivepoint（0.025,80,40,10000,1e-11）;
k=3952;
figure
surf（x,y,u）;xlswrite（'D:
\define.xls',e,'data1'）;%保存数据
xlabel（'x'）;ylabel（'y'）;zlabel（'u'）;
title（'五点差分法解椭圆型偏微分方程'）;
figure
surf（x,y,e）;xlswrite（'D:
\wucha.xls',e,'data2'）;%保存数据
title（'步长为0.025时的误差'）;
figure
surf（x,y,p）;xlswrite（'D:
\jinshi.xls',e,'data32'）;%保存数据
title（'精确解曲面图'）;
结果输出：
精确解
近似解
误差
2.02E-05
2.02E-05
2.02E-05
4.04E-05
4.04E-05
4.04E-05
6.02E-05
6.02E-05
6.02E-05
7.97E-05
7.97E-05
7.97E-05
9.87E-05
9.87E-05
9.87E-05
0.000117
0.000117
0.000117
0.000135
0.000135
0.000135
0.000152
0.000152
0.000152
0.000168
0.000168
0.000168
0.000182
0.000182
0.000182
0.000196
0.000196
0.000196
0.000209
0.000209
0.000209
0.00022
0.00022
0.00022
0.00023
0.00023
0.00023
0.000238
0.000238
0.000238
0.000245
0.000245
0.000245
0.000251
0.000251
0.000251
0.000255
0.000255
0.000255
0.000257
0.000257
0.000257
0.000258
0.000258
0.000258
0.000257
0.000257
0.000257
0.000255
0.000255
0.000255
升级会员 