大学本科常微分实验报告.docx
《大学本科常微分实验报告.docx》由会员分享，可在线阅读，更多相关《大学本科常微分实验报告.docx（25页珍藏版）》请在冰豆网上搜索。
大学本科常微分实验报告
福建农林大学计算机与信息学院
（数学类课程）
课程实习报告
课程名称：
常微分方程课程实习
实习题目：
常微分方程数值求解问题的实习
姓名：
系：
信息与计算科学
专业：
信息与计算科学
年级：
2010
学号：
指导教师：
职称：
讲师
2011年12月1日
福建农林大学计算机与信息学院数学类
课程实习报告结果评定
评语：
成绩：
指导教师签字：
评定日期：
目录
1.实习的目的和任务1
2.实习要求1
3.实习地点1
4.主要仪器设备1
5.实习内容1
5.1用不同格式对同一个初值问题的数值求解及其分析……………………..1
5.1.1求精确解1
5.1.2用欧拉法求解3
5.1.3用改进欧拉法求解5
5.1.4用4级4阶龙格—库塔法求解7
5.1.5问题讨论与分析………………………………………………………………………10
5.2一个算法不同不长求解同一个初值问题及其分析…………………………………..13
5.3洛伦茨方程模拟混沌现象……………………………………………………………18
6.结束语21
参考文献21
常微分方程课程实习
1.实习的目的和任务
目的：
通过课程实习能够应用MATLAB软来计算微分方程（组）的数值解；了解常微分方程数值解。
任务：
通过具体的问题，利用MATLAB软件来计算问题的结果，分析问题的结论。
2.实习要求
能够从案例的自然语言描述中，抽象出其中的数学模型；能够熟练应用所学的数值解计算方法；能够熟练使用MATLAB软件；对常微分方程数值解有所认识，包括对不同算法有所认识和对步长有所认识。
3.实习地点
南2#425
4.主要仪器设备
计算机
MicrosoftWindows7
MatlabR2009a
5.实习内容
5.1用欧拉方法，改进欧拉方法，4阶龙格—库塔方法分别求下面微分方程的初值dy/dx=y*cos（2*x）y（0）=1x∈[0，2]
5.1.1求精确解
变量分离方程情形:
形如
的方程,这里
分别是
的连续函数.如果
我们可将方程改写成
这样,变量就”分离”开来了,两边同时积分即可:
为任意常数.
用变量分离法可求出其精确为：
y=exp（0.5*sin（2*x））
5.1.1程序代码：
>>x=0:
0.1:
2;
>>y=exp（0.5*sin（2*x））
>>plot（x,y,'rs-'）;
>>Data=[x',y']
结果及图像：
y=
Columns1through6
1.00001.10441.21501.32621.43141.5231
Columns7through12
1.59361.63681.64841.62731.57561.4982
Columns13through18
1.40181.29401.18231.07310.97120.8801
Columns19through21
0.80150.73640.6850
Data=
01.0000
0.10001.1044
0.20001.2150
0.30001.3262
0.40001.4314
0.50001.5231
0.60001.5936
0.70001.6368
0.80001.6484
0.90001.6273
1.00001.5756
1.10001.4982
1.20001.4018
1.30001.2940
1.40001.1823
1.50001.0731
1.60000.9712
1.70000.8801
1.80000.8015
1.90000.7364
2.00000.6850
5.1.2用欧拉法求解
设常微分方程的初始问题
有唯一解。
则由欧拉法求初值问题
（1），
（2）的数值解的计算公式为：
（
）
程序如下：
建立函数文件cwf1.m
function[x,y]=cwf1（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（'y*cos（2*x）'）;
>>[x,y]=cwf1（fun,[0,2],1,0.1）;
>>[x,y]
>>plot（x,y,'r*-'）
结果及图像：
ans=
01.0000
0.10001.1000
0.20001.2078
0.30001.3191
0.40001.4279
0.50001.5274
0.60001.6099
0.70001.6683
0.80001.6966
0.90001.6917
1.00001.6532
1.10001.5844
1.20001.4912
1.30001.3812
1.40001.2629
1.50001.1439
1.60001.0306
1.70000.9278
1.80000.8381
1.90000.7629
2.00000.7026
5.1.3用改进欧拉法求解：
计算公式为：
当然也可以迭代多次：
程序如下：
建立函数文件cwf2.m
function[x,y]=cwf2（fun,x_span,y0,h）
x=x_span
（1）:
h:
x_span
（2）;
y
（1）=y0;
forn=1:
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;
end
x=x';y=y';
在MATLAB输入以下程序：
>>clearall
>>fun=inline（'y*cos（2*x）'）;
>>[x,y]=cwf2（fun,[0,2],1,0.1）;
>>[x,y]
>>plot（x,y,'b+-'）
结果及图像：
ans=
01.0000
0.10001.1039
0.20001.2138
0.30001.3244
0.40001.4290
0.50001.5201
0.60001.5902
0.70001.6330
0.80001.6445
0.90001.6234
1.00001.5720
1.10001.4949
1.20001.3991
1.30001.2920
1.40001.1810
1.50001.0724
1.60000.9711
1.70000.8803
1.80000.8021
1.90000.7373
2.00000.6859
5.1.4用4阶龙格—库塔求解
四级四阶的龙格—库塔的一般算式
公式的截断误差阶为
。
经典四级四阶龙格—库塔格式取定
，则得
这是最为著名的经典四级四阶龙格—库塔格式。
程序如下：
建立函数文件cwf3.m
function[x,y]=cwf3（fun,x_span,y0,h）
x=x_span
（1）:
h:
x_span
（2）;
y
（1）=y0;
forn=1:
length（x）-1
k1=feval（fun,x（n）,y（n））;
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;
end
x=x';y=y';
在MATLAB输入以下程序：
>>clearall;
>>fun=inline（'y*cos（2*x）'）;
>>[x,y]=cwf3（fun,[0,2],1,0.1）;
>>[x,y]
>>plot（x,y,'b+-'）
结果及其图象：
ans=
01.0000
0.10001.1044
0.20001.2150
0.30001.3262
0.40001.4314
0.50001.5231
0.60001.5936
0.70001.6368
0.80001.6484
0.90001.6273
1.00001.5756
1.10001.4982
1.20001.4018
1.30001.2940
1.40001.1823
1.50001.0731
1.60000.9712
1.70000.8801
1.80000.8015
1.90000.7364
2.00000.6850
5.1.5问题讨论与分析
由以上数值分析结果绘制表格：
精确解
欧拉方法
改进的欧拉方法
四阶龙格-库塔方法
xi
yi
yi
误差
yi
误差
yi
误差
0
1
1
0
1
0
1
0
0.1
1.1044
1.1
0.0044
1.1039
0.0005
1.1044
0
0.2
1.215
1.2078
0.0072
1.2138
0.0012
1.215
0
0.3
1.3262
1.3191
0.0071
1.3244
0.0018
1.3262
0
0.4
1.4314
1.4279
0.0035
1.429
0.0024
1.4314
0
0.5
1.5231
1.5274
-0.0043
1.5201
0.003
1.5231
0
0.6
1.5936
1.6099
-0.0163
1.5902
0.0034
1.5936
0
0.7
1.6368
1.6683
-0.0315
1.633
0.0038
1.6368
0
0.8
1.6484
1.6966
-0.0482
1.6445
0.0039
1.6484
0
0.9
1.6273
1.6917
-0.0644
1.6234
0.0039
1.6273
0
1.0
1.5756
1.6532
-0.0776
1.572
0.0036
1.5756
0
1.1
1.4982
1.5844
-0.0862
1.4949
0.0033
1.4982
0
1.2
1.4018
1.4912
-0.0894
1.3991
0.0027
1.4018
0
1.3
1.294
1.3812
-0.0872
1.292
0.002
1.294
0
1.4
1.1823
1.2629
-0.0806
1.181
0.0013
1.1823
0
1.5
1.0731
1.1439
-0.0708
1.0724
0.0007
1.0731
0
1.6
0.9712
1.0306
-0.0594
0.9711
0.0001
0.9712
0
1.7
0.8801
0.9278
-0.0477
0.8803
-0.0002
0.8801
0
1.8
0.8015
0.8381
-0.0366
0.8021
-0.0006
0.8015
0
1.9
0.7364
0.7629
-0.0265
0.7373
-0.0009
0.7364
0
2
0.685
0.7026
-0.0176
0.6859
-0.0009
0.685
0
在MATLAB输入以下程序：
>>x=0:
0.1:
2;
>>y1=[1.00001.10441.21501.32621.43141.52311.59361.63681.64841.62731.57561.49821.40181.29401.18231.07310.97120.88010.80150.73640.6850];
>>y2=[1.00001.10001.20781.31911.42791.52741.60991.66831.69661.69171.65321.58441.49121.38121.26291.14391.03060.92780.83810.76290.7026];
>>y3=[1.00001.10391.21381.32441.42901.52011.59021.63301.64451.62341.57201.49491.39911.29201.18101.07240.97110.88030.80210.73730.6859];
>>y4=[1.00001.10441.21501.32621.43141.52311.59361.63681.64841.62731.57561.49821.40181.29401.18231.07310.97120.88010.80150.73640.6850];
>>plot（x,y1,'r+-'）
>>holdon,plot（x,y2,'b-'）
>>plot（x,y1,'r+-'）
>>holdon,plot（x,y3,'b-'）
>>plot（x,y1,'r+-'）
>>holdon,plot（x,y4,'b-'）
结果如下：
精确解与欧拉方法的比较
精确解与改进欧拉方法的比较
精确解与4阶龙格—库塔方法比较
结果分析：
由上表和图可以看出欧拉法误差最大，而改进欧拉和龙格—库塔方法误差相对较小，并且龙格—库塔方法误差最小且大部分值都跟精确值相同。
由欧拉图与精确图相比可清晰看到，随着x的增加，函数值与精确值的偏差越来越大。
5.2选择用欧拉方法取不同步长分别求下面微分方程的初值dy/dx=y*cos（2*x）
y（0）=1，x∈[0，2]。
程序如下：
建立函数文件cwf1.m
function[x,y]=cwf1（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';
①当h=0.05时
在MATLAB输入以下程序：
>>clearall
>>fun=inline（'y*cos（2*x）'）;
>>[x,y]=cwf1（fun,[0,2],1,0.05）;
>>[x,y]
>>plot（x,y,'r*-'）
结果及图像：
ans=
01.0000
0.05001.0500
0.10001.1022
0.15001.1563
0.20001.2115
0.25001.2673
0.30001.3229
0.35001.3775
0.40001.4301
0.45001.4800
0.50001.5260
0.55001.5672
0.60001.6027
0.65001.6318
0.70001.6536
0.75001.6677
0.80001.6735
0.85001.6711
0.90001.6603
0.95001.6415
1.00001.6149
1.05001.5813
1.10001.5414
1.15001.4961
1.20001.4462
1.25001.3929
1.30001.3371
1.35001.2798
1.40001.2220
1.45001.1644
1.50001.1079
1.55001.0530
1.60001.0004
1.65000.9505
1.70000.9036
1.75000.8599
1.80000.8196
1.85000.7829
1.90000.7497
1.95000.7200
2.00000.6939
②当h=0.1时
在MATLAB输入以下程序：
>>clearall
>>fun=inline（'y*cos（2*x）'）;
>>[x,y]=cwf1（fun,[0,2],1,0.1）;
>>[x,y]
>>plot（x,y,'r*-'）
结果及图像：
ans=
01.0000
0.10001.1000
0.20001.2078
0.30001.3191
0.40001.4279
0.50001.5274
0.60001.6099
0.70001.6683
0.80001.6966
0.90001.6917
1.00001.6532
1.10001.5844
1.20001.4912
1.30001.3812
1.40001.2629
1.50001.1439
1.60001.0306
1.70000.9278
1.80000.8381
1.90000.7629
2.00000.7026
③当h=0.5时
在MATLAB输入以下程序：
>>clearall
>>fun=inline（'y*cos（2*x）'）;
>>[x,y]=cwf1（fun,[0,2],1,0.5）;
>>[x,y]
>>plot（x,y,'r*-'）
结果及图像：
ans=
01.0000
0.50001.5000
1.00001.9052
1.50001.5088
2.00000.7619
结果分析：
根据以上的结果和图像与精确解的比较可知，步长越小，用欧拉方法所求的值就越接近精确解。
5.3数值求解Lorenz方程，模拟混沌现象
洛伦茨方程如下：
将x，y，z表示y
（1），y
（2），y（3），即列向量y中三个分量。
程序如下：
（1）建立自定义函数
用edit命令建立自定义函数名为lorenz.m的内容为：
functiondy=lorenz（t,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/3*y（3）
（2）用edit命令建立一个命令文件xian.m的内容为
[t,y]=ode45（'lorenz',[0，30],[12,2,9]）;
%表示在0~30秒内求解，在零时刻设y
（1）为12，y
（2）为2，y（3）为9.
plot（t,y（:
1））%显示y
（1）即x与时间的关系图
pause%暂停
plot（t,y（:
2））%显示y
（2）即y与时间的关系图
pause
plot（t,y（:
3））%显示y（3）即z与时间的关系图
pause
plot3（y（:
1）,y（:
2）,y（:
3））;%显示三分量的关系图
view（[20,42]）;%设定一个较好的观察角度
（3）求解结果
在matlab窗口中执行xian.m文件，数据偏多就不显示了，其图像如下：
y
（1）即x与时间的关系图
y
（2）即y与时间的关系图
y（3）即z与时间的关系图
显示三分量的关系图
（4）用edit命令建立一个命令文件wea.m的内容为
[t,y]=ode45（'lorenz',[0，30],[3,7,18]）;
%表示在0~30秒内求解，在零时刻设y
（1）为12，y
（2）为2，y（3）为9.
plot（t,y（:
1））%显示y
（1）即x与时间的关系图
pause%暂停
plot（t,y（:
2））%显示y
（2）即y与时间的关系图
pause
plot（t,y（:
3））%显示y（3）即z与时间的关系图
pause
plot3（y（:
1）,y（:
2）,y（:
3））;%显示三分量的关系图
view（[20,42]）;%设定一个较好的观察角度
在matlab窗口中执行wea.m文件，其图像如下：
显示三分量的关系图
结果分析：
当参数特定选定时，初值稍有变化，函数图像的绕法就会发生变化，也就是说微分系统对初值的敏感度相当的大，这就形成了所谓的“混沌”。
6、结束语
通过这段期间的实习，我基本掌握了MATLAB语言的一些基本操作，特别是绘图功能及编程能力.由于有些知识以前没接触过，通过实习提高了我的临时自学能力，同时提高了自己抽象思维的能力和从数学模型中理解问题所在。
在这次的实习中，让我明白了网络的方便与简洁，对于不懂的问题可以及时得到解决，同时我也知道了学校的图书馆原来有那么多的学习资源，以后要多去图书馆吸收知识，充实自己。
参考文献：
[1]张圣勤.《MATLAB7.0实用教程》.北京:
机械工业出版社，2004
[2]姜健飞、胡良剑、唐俭.《数值分析及其MATLAB实验》．上海：
科学出版社，2004
[3]许波，刘征.MATLAB工程数学应用［M］.北京：
清华大学出版社，2000
[4]刘华杰.分形艺术（电子版［M］.长沙：
湖南电子音像出版社，1997
升级会员 