福建农林大学常微分课程论文Word格式文档下载.docx
《福建农林大学常微分课程论文Word格式文档下载.docx》由会员分享，可在线阅读，更多相关《福建农林大学常微分课程论文Word格式文档下载.docx（23页珍藏版）》请在冰豆网上搜索。
6.结束语19
参考文献19
1.实习的目的和任务
目的：
通过课程实习能够应用MATLAB软来计算微分方程（组）的数值解；
了解常微分方程数值解。
任务：
通过具体的问题，利用MATLAB软件来计算问题的结果，分析问题的结论。
2.实习要求
能够从案例的自然语言描述中，抽象出其中的数学模型；
能够熟练应用所学的数值解计算方法；
能够熟练使用MATLAB软件；
对常微分方程数值解有所认识，包括对不同算法有所认识和对步长有所认识。
3.实习地点
学生宿舍
4.主要仪器设备
计算机、
MicrosoftWindowsXP
Matlab7.6
5.实习内容
5.1数值求解Lorenz方程，模拟混沌现象
Lorenz方程：
，其中a=10,b=8/3,c=28.
5.1.1数值求解方程
程序代码：
建立函数文件
functionxdot=lorenzeq（t,x）
a=10;
r=28;
b=8/3;
xdot=[a*（x
（2）-x
（1））;
r*x
（1）-x
（2）-x
（1）*x
（2）;
x
（1）*x
（2）-b*x（3）];
调用函数
globalabc
c=28;
%定义参数
lorenz=@（t,x）[a*（x
（2）-x
（1））;
c*x
（1）-x
（1）*x（3）-x
（2）;
x
（1）*x
（2）-b*x（3）];
%定义函数
[T,X]=ode45（lorenz,[0,20],[0;
1;
2]）;
%数值法解微分方程
Data=[T,X]
plot3（X（:
1）,X（:
2）,X（:
3）,'
m'
）%绘图
view（-20,60）;
%设置视角
xlabel（'
x'
）;
ylabel（'
y'
zlabel（'
z'
%标记坐标轴
Data=
001.00002.0000
0.00010.00120.99991.9994
0.00020.00250.99981.9987
0.00040.00370.99961.9980
0.00050.00500.99951.9974
...................................
19.9735-10.00070.259438.1095
19.9801-9.32150.881337.4048
19.9867-8.64891.413036.6808
19.9934-7.98811.861135.9479
20.0000-7.34352.233135.2142
运行则得到参数a=10，b=8/3，c=28自变量t取值为[0,20]的Lorenz方程图像（图1）和数据，图中出现一个吸引子，整个图像由螺线型轨道构成。
图1
5.1.2数值解对初值的敏感性
将初值x（0）=0稍微改变下，取x（0）=0.001和x（0）=0进行比较
程序代码：
globalabc
a=10;
c*x
（1）-x
（1）*x（3）-x
（2）;
%定义函数
subplot（2,2,1）;
[T,X]=ode45（lorenz,[0,100],[0;
plot（T,X（:
1）,'
m-'
xlabel（'
t'
holdon
[T,X]=ode45（lorenz,[0,100],[0.001;
b'
）
holdoff
subplot（2,2,2）;
2）,'
subplot（2,2,3）;
图2
由图2可知：
当初始条件x（0）=0稍微改变为x（0）=0.001其他两个初始条件不变时，Lorenz方程的三个自变量的数值解变化很大，同理可得另外两个初始条件稍微变化后的结果，这里不详细讨论，从而可以简单说明Lorenz方程的数值解对初值很敏感。
5.1.3模拟混沌现象
将t取不同值分别绘出不同值的函数图像
:
lorenz=@（t,x）[a*（x
（2）-x
（1））;
[T,X]=ode45（lorenz,[0,10],[2;
3;
4]）;
holdon
subplot（2,2,2）;
[T,X]=ode45（lorenz,[0,20],[2;
subplot（2,2,3）;
[T,X]=ode45（lorenz,[0,100],[2;
subplot（2,2,4）;
[T,X]=ode45（lorenz,[0,500],[2;
holdoff
图3
由图3可知：
方程有两个平衡点，当初值和参数不变时，随着计算次数的增加，相空间曲线既不趋近于平衡点1，也不趋近于平衡点2，始终在两个平衡点一定范围内无限反复旋转，此时方程进入混沌状态，不管t值取多大，方程依然处于混沌状态，根本无法预测绕平衡点的旋转次数，从而最终得出Lorenz方程对初值非常敏感。
5.2用欧拉方法，改进欧拉方法，4阶龙格—库塔方法分别求下面微分方程的初值dy/dx=（1+y^2）*sin（x）y（0）=1x∈[-1.2，1.2]
5.2.1求精确解
首先可以求得其精确解为：
y=cot（cos（x）-1+1/4*pi）
由matlab可得函数y的图像（图4），由图知该函数在区间上有间断点，但可以取其中一个连续区间[-1.2,1.2]进行研究。
图4
>
x=-1.2:
0.1:
1.2;
y=cot（cos（x）-1+1/4*pi）;
plot（x,y,'
b*-'
Data=[x'
y'
]
Data=
-1.20006.7186
-1.10004.1042
-1.00002.9610
-0.90002.3198
-0.80001.9110
-0.70001.6302
-0.60001.4285
-0.50001.2806
-0.40001.1718
-0.30001.0936
-0.20001.0407
-0.10001.0100
01.0000
0.10001.0100
0.20001.0407
0.30001.0936
0.40001.1718
0.50001.2806
0.60001.4285
0.70001.6302
0.80001.9110
0.90002.3198
1.00002.9610
1.10004.1042
1.20006.7186
图5
5.2.2用欧拉法求解
程序如下：
建立函数文件cwfa1.m
function[x,y]=cwfa1（fun,x_span,y0,h）
x=x_span
（1）:
h:
x_span
（2）;
y
（1）=y0;
forn=1:
length（x）-1
y（n+1）=y（n）+h*feval（fun,x（n）,y（n））;
end
x=x'
;
y=y'
在MATLAB输入以下程序：
clearall
fun=inline（'
（1+y^2）*sin（x）'
[x,y]=cwfa1（fun,[-1.2,1.2],7,0.1）;
[x,y]
r*-'
）
图6
5.2.3用改进欧拉法求解：
建立函数文件cwfa2.m
function[x,y]=cwfa2（fun,x_span,y0,h）
length（x）-1
k1=feval（fun,x（n）,y（n））;
y（n+1）=y（n）+h*k1;
k2=feval（fun,x（n+1）,y（n+1））;
y（n+1）=y（n）+h*（k1+k2）/2;
（1+y^2）*sin（x）'
[x,y]=cwfa2（fun,[-1.2,1.2],7,0.1）;
k+-'
图7
5.2.4用4阶龙格—库塔求解
建立函数文件cwfa3.m
function[x,y]=cwfa3（fun,x_span,y0,h）
k2=feval（fun,x（n）+h/2,y（n）+h/2*k1）;
k3=feval（fun,x（n）+h/2,y（n）+h/2*k2）;
k4=feval（fun,x（n+1）,y（n）+h*k3）;
y（n+1）=y（n）+h*（k1+2*k2+2*k3+k4）/6;
clearall;
（1+y^2）*sin（x）'
[x,y]=cwfa3（fun,[-1.2,1.2],7,0.1）;
[x,y]
plot（x,y,'
g+-'
图8
5.2.5问题讨论与分析
由以上数值分析结果绘制表格：
精确解
欧拉方法
改进的欧拉方法
四阶龙格-库塔方法
xi
yi
误差
-1.2
6.7186
7
-0.2814
-1.1
4.1042
2.3398
1.7644
4.3814
-0.2772
4.2039
-0.0997
-1
2.961
1.7628
1.1982
3.159
-0.198
3.015
-0.054
-0.9
2.3198
1.4172
0.9026
2.4622
-0.1424
2.3549
-0.0351
-0.8
1.911
1.1815
0.7295
2.019
-0.108
1.9365
-0.0255
-0.7
1.6302
1.0096
0.6206
1.7165
-0.0863
1.6502
-0.02
-0.6
1.4285
0.8795
0.549
1.5008
-0.0723
1.4452
-0.0167
-0.5
1.2806
0.7794
0.5012
1.3434
-0.0628
1.295
-0.0144
-0.4
1.1718
0.7023
0.4695
1.2283
-0.0565
1.1848
-0.013
-0.3
1.0936
0.6442
0.4494
1.1457
-0.0521
1.1056
-0.012
-0.2
1.0407
0.6024
0.4383
1.0901
-0.0494
1.0521
-0.0114
-0.1
1.01
0.5753
0.4347
1.0579
-0.0479
1.0211
-0.0111
1
0.562
0.438
1.0473
-0.0473
1.0109
-0.0109
0.1
0.448
1.0578
-0.0478
0.2
0.5751
0.4656
1.0899
-0.0492
0.3
0.6016
0.492
1.1454
-0.0518
0.4
0.6418
0.53
1.2277
-0.0559
0.5
0.6968
0.5838
1.3426
-0.062
0.6
0.768
0.6605
1.4996
-0.0711
0.7
0.8578
0.7724
1.7147
-0.0845
0.8
0.9696
0.9414
2.0165
-0.1055
0.9
1.1088
1.211
2.4592
-0.1394
2.3548
-0.035
1.2834
1.6776
3.1589
-0.1979
3.0148
-0.0538
1.1
1.5062
2.598
4.4081
-0.3039
4.2027
-0.0985
1.2
1.7975
4.9211
7.1733
-0.4547
6.9689
-0.2503
将三种方法图像与精确解作比较
y=cot（cos（x）-1+1/4*pi）;
plot（x,y,'
fun=inline（'
[x,y]=cwfa2（fun,[-1.2,1.2],7,0.1）;
[x,y]=cwfa1（fun,[-1.2,1.2],7,0.1）;
rx-'
[x,y]=cwfa3（fun,[-1.2,1.2],7,0.1）;
plot（x,y,'
g.-'
图9
分析和结论
由上表和图可以看出欧拉法误差最大，而改进欧拉和龙格—库塔方法误差相对较小，并且龙格—库塔方法误差最小且几乎都跟精确值相同。
由欧拉图与精确图相比可清晰看到，随着x由0向两边增加，函数值与精确值的偏差越来越大。
5.3用改进欧拉方法取不同步长分别求下面微分方程的初值dy/dx=（1+y^2）*sin（x），y（0）=1x∈[-1.2，1.2]，并对结果进行分析说明，给出你的结论。
5.3.1取步长
h分别取0.12、0.14、0.16、0.18、0.2
c-'
holdon;
[x,y]=cwfa2（fun,[-1.2,1.2],7,0.12）;
[x,y]
g*-'
[x,y]=cwfa2（fun,[-1.2,1.2],7,0.14）;
r+:
'
[x,y]=cwfa2（fun,[-1.2,1.2],7,0.16）;
[x,y]=cwfa2（fun,[-1.2,1.2],7,0.18）;
k+:
[x,y]=cwfa2（fun,[-1.2,1.2],7,0.2）;
m*-'
图10
图11
图12
图13
5.3.2分析和结论
由于原本比较图像的线比较密集，采取分别对图像的左部（图13）、中部（图11）、右部（图12）进行放大，由放大的图像可知h=0.12的精确度最高，最差是h=0.2，随着h的增加所取点数的减少，所得的图像与精确解的相差越来越大；
h取得越小精确度越高。
6.结束语
在没做这份实习报告之前对Lorenz方程和数值解只是理论上的认识。
Lorenz方程没有解析解，只有数值解。
Lorenz方程还有对称性、z轴是不变集、耗散性和吸引性的性质。
此实习只研究方程对初值的敏感性和混沌现象，自己在网上下载了几篇关于这方面的论文和实验报告，借鉴了其中一些好的方法，加上自己的想法整个实习形成了自己的一套思路，操作的过程中要用到几个关键程序代码，自己通过在网上找资料和查阅书籍解决了编程问题，图的绘制和修饰自己在平时的操作练习中也已初步掌握。
从课外书籍了解到：
对于Lorenz方程，不同参数的设置可以得到不同的状态。
对欧拉法、改进欧拉法和四阶龙格-库塔方法，从实习的结果认识到对于一般的微分方程用欧拉法作数值求解就可以达到足够精确了，前提是步长要取得够小，一般数值求解用的改进欧拉法。
从寻找方程的过程中认识到这三种并不能适应于所有一阶微分方程，就连四阶龙格-库塔方法也一样。
这次实习课程设计很大程度上提高了我的动手和实习能力，让我认识到了实践和理论相结合的重要性。
参考文献
[1]刘卫国.《MATLAB程序设计与应用》.北京:
高等教育出版社，2006
[2]王高雄、周之铭、朱思铭、王寿松.《常微分方程》．北京:
[3]黄振侃.《数值计算-常微分方程数值解》．北京:
北京工业大学，2006
[4]（美）赫希（Hirsch，M.W）、（美）斯梅尔.《（全新正版）微分方程动力系统和混沌导论第二版》．北京:
世界图书出版社，2007
升级会员 