常微分方程课程实习.docx
《常微分方程课程实习.docx》由会员分享，可在线阅读，更多相关《常微分方程课程实习.docx（28页珍藏版）》请在冰点文库上搜索。
常微分方程课程实习
福建农林大学计算机与信息学院
（数学类课程）
课程实习报告
课程名称：
常微分方程课程实习
实习题目：
常微分方程数值求解问题的实习
姓名：
系：
应用数学
专业：
数学与应用数学
年级：
2008
学号：
指导教师：
职称：
讲师
2009年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问题讨论与分析………………………………………………………………………9
5.2一个算法不同不长求解同一个初值问题及其分析…………………………………..13
6.结束语22
参考文献22
常微分方程课程实习
1.实习的目的和任务
目的：
通过课程实习能够应用MATLAB软来计算微分方程（组）的数值解；了解常微分方程数值解。
任务：
通过具体的问题，利用MATLAB软件来计算问题的结果，分析问题的结论。
2.实习要求
能够从案例的自然语言描述中，抽象出其中的数学模型；能够熟练应用所学的数值解计算方法；能够熟练使用MATLAB软件；对常微分方程数值解有所认识，包括对不同算法有所认识和对步长有所认识。
3.实习地点
数学实验室
4.主要仪器设备
计算机、
MicrosoftWindows7
MatlabR2009a
5.实习内容
5.1用欧拉方法，改进欧拉方法，4阶龙格—库塔方法分别求下面微分方程的初值dy/dx=y-x^2+1y（0）=0.5x∈[0，2]
5.1.1求精确解
首先可以求得其精确解为：
y=2*x-exp（x）/2+x^2+1
5.1.1程序代码：
>>x=0:
0.1:
2;
y=2*x-exp（x）/2+x.^2+1
plot（x,y,'r*-'）
Data=[x',y']
y=
Columns1through17
0.50000.65740.82931.01511.21411.42561.64891.88312.12722.38022.64092.90793.17993.45543.73244.00924.2835
Columns18through21
4.55304.81525.06715.3055
Data=
00.5000
0.10000.6574
0.20000.8293
0.30001.0151
0.40001.2141
0.50001.4256
0.60001.6489
0.70001.8831
0.80002.1272
0.90002.3802
1.00002.6409
1.10002.9079
1.20003.1799
1.30003.4554
1.40003.7324
1.50004.0092
1.60004.2835
1.70004.5530
1.80004.8152
1.90005.0671
2.00005.3055
5.1.2用欧拉法求解
程序如下：
建立函数文件Euler.m
function[x,y]=Euler（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-x^2+1'）;
>>[x,y]=Euler（fun,[0,2],0.5,0.1）;
>>[x,y]
>>plot（x,y,'b*-'）
结果及其图象：
ans=
00.5000
0.10000.6500
0.20000.8140
0.30000.9914
0.40001.1815
0.50001.3837
0.60001.5971
0.70001.8208
0.80002.0538
0.90002.2952
1.00002.5438
1.10002.7981
1.20003.0569
1.30003.3186
1.40003.5815
1.50003.8437
1.60004.1030
1.70004.3573
1.80004.6040
1.90004.8405
2.00005.0635
5.1.3用改进欧拉法求解：
程序如下：
建立函数文件gjEuler.m
function[x,y]=gjEuler（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-x^2+1'）;
>>[x,y]=gjEuler（fun,[0,2],0.5,0.1）;
>>[x,y]
>>plot（x,y,'m+-'）
结果及其图象：
00.5000
0.10000.6570
0.20000.8284
0.30001.0137
0.40001.2122
0.50001.4232
0.60001.6459
0.70001.8794
0.80002.1228
0.90002.3750
1.00002.6348
1.10002.9010
1.20003.1720
1.30003.4464
1.40003.7223
1.50003.9978
1.60004.2708
1.70004.5390
1.80004.7996
1.90005.0499
2.00005.2866
5.1.4用4阶龙格—库塔求解
程序如下：
建立函数文件Runge_Kutta.m
function[x,y]=Runge_Kutta（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-x^2+1'）;
>>[x,y]=Runge_Kutta（fun,[0,2],0.5,0.1）;
>>[x,y]
>>plot（x,y,'g+-'）
结果及其图象：
ans=
ans=
00.5000
0.10000.6574
0.20000.8293
0.30001.0151
0.40001.2141
0.50001.4256
0.60001.6489
0.70001.8831
0.80002.1272
0.90002.3802
1.00002.6409
1.10002.9079
1.20003.1799
1.30003.4553
1.40003.7324
1.50004.0092
1.60004.2835
1.70004.5530
1.80004.8152
1.90005.0670
2.00005.3055
5.1.5问题讨论与分析
由以上数值分析结果绘制表格：
精确解
欧拉方法
改进的欧拉方法
四阶龙格-库塔方法
xi
yi
yi
误差
yi
误差
yi
误差
0
0.5000
0.5000
0
0.5000
0
0.5000
0
0.1
0.6574
0.6500
-0.0074
0.6570
-0.0004
0.6574
0
0.2
0.8293
0.8140
-0.0153
0.8284
-0.0009
0.8293
0
0.3
1.0151
0.9914
-0.0237
1.0137
-0.0014
1.0151
0
0.4
1.2141
1.1815
-0.0326
1.2122
-0.0019
1.2141
0
0.5
1.4256
1.3837
-0.0419
1.4232
-0.0024
1.4256
0
0.6
1.6489
1.5971
-0.0518
1.6459
-0.0030
1.6489
0
0.7
1.8831
1.8208
-0.0623
1.8794
-0.0037
1.8831
0
0.8
2.1272
2.0538
-0.0734
2.1228
-0.0044
2.1272
0
0.9
2.3802
2.2952
-0.0850
2.3750
-0.0052
2.3802
0
1.0
2.6409
2.5438
-0.0971
2.6348
-0.0061
2.6409
0
1.1
2.9079
2.7981
-0.1098
2.9010
-0.0069
2.9079
0
1.2
3.1799
3.0569
-0.1230
3.1720
-0.0079
3.1799
0
1.3
3.4554
3.3186
-0.1368
3.4464
-0.0090
3.4553
-0.0001
1.4
3.7324
3.5815
-0.1509
3.7223
-0.0101
3.7324
0
1.5
4.0092
3.8437
-0.1655
3.9978
-0.0114
4.0092
0
1.6
4.2835
4.1030
-0.1805
4.2708
-0.0127
4.2835
0
1.7
4.5530
4.3573
-0.1957
4.5390
-0.0140
4.5530
0
1.8
4.8152
4.6040
-0.2112
4.7996
-0.0156
4.8152
0
1.9
5.0671
4.8405
-0.2266
5.0499
-0.0172
5.0670
-0.0001
2.0
5.3055
5.0635
-0.2420
5.2866
-0.0189
5.3055
0
x=0:
0.1:
2;
y1=[0.50000.65740.82931.01511.21411.42561.64891.88312.12722.38022.64092.90793.17993.45543.73244.00924.28354.55304.81525.06715.3055];
>>y1=[0.50000.65740.82931.01511.21411.42561.64891.88312.12722.38022.64092.90793.17993.45543.73244.00924.28354.55304.81525.06715.3055];
>>y2=[0.50000.65000.81400.99141.18151.38371.59711.82082.05382.29522.54382.79813.05693.31863.58153.84374.10304.35734.60404.8405.0635];
>>y3=[0.50000.65700.82841.01371.21221.42321.64591.87942.12282.37502.63482.90103.17203.44643.72233.99784.27084.53904.79965.04995.2866];
>>y4=[0.50000.65740.82931.01511.21411.42561.64891.88312.12722.38022.64092.90793.17993.45533.73244.00924.28354.55304.81525.06705.3055];
>>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-'）
由上表和图可以看出欧拉法误差最大，而改进欧拉和龙格—库塔方法误差相对较小，并且龙格—库塔方法误差最小且大部分值都跟精确值相同。
由欧拉图与精确图相比可清晰看到，随着x的增加，函数值与精确值的偏差越来越大。
5.2选择用欧拉方法，改进欧拉方法，4阶龙格—库塔方法之一取不同步长分别求下面微分方程的初值dy/dx=y-x^2+1y（0）=0.5x∈[0，2]，并对结果进行分析说明，给出你的结论。
用改进欧拉方法取不同步长分别求微分方程的初值dy/dx=y-x^2+1y（0）=0.5x∈[0，2]。
一．步长0.05
在MATLAB输入以下程序：
>>clearall
>>fun=inline（'y-x^2+1'）;
>>[x,y]=gjEuler（fun,[0,2],0.5,0.05）;
>>[x,y]
>>plot（x,y,'r*-'）
结果及其图像：
ans=
00.5000
0.05000.5768
0.10000.6573
0.15000.7414
0.20000.8291
0.25000.9202
0.30001.0147
0.35001.1126
0.40001.2136
0.45001.3178
0.50001.4250
0.55001.5352
0.60001.6482
0.65001.7639
0.70001.8822
0.75002.0030
0.80002.1261
0.85002.2514
0.90002.3789
0.95002.5082
1.00002.6393
1.05002.7720
1.10002.9061
1.15003.0415
1.20003.1779
1.25003.3152
1.30003.4531
1.35003.5913
1.40003.7298
1.45003.8682
1.50004.0063
1.55004.1437
1.60004.2803
1.65004.4156
1.70004.5494
1.75004.6814
1.80004.8112
1.85004.9384
1.90005.0627
1.95005.1836
2.00005.3007
二．步长0.1
在MATLAB输入以下程序：
>>clearall
>>fun=inline（'y-x^2+1'）;
>>[x,y]=gjEuler（fun,[0,2],0.5,0.1）;
>>[x,y]
>>plot（x,y,'m+-'）
结果及其图象：
00.5000
0.10000.6570
0.20000.8284
0.30001.0137
0.40001.2122
0.50001.4232
0.60001.6459
0.70001.8794
0.80002.1228
0.90002.3750
1.00002.6348
1.10002.9010
1.20003.1720
1.30003.4464
1.40003.7223
1.50003.9978
1.60004.2708
1.70004.5390
1.80004.7996
1.90005.0499
2.00005.2866
三．步长0.15
在MATLAB输入以下程序：
>>clearall
>>fun=inline（'y-x^2+1'）;
>>[x,y]=gjEuler（fun,[0,2],0.5,0.15）;
>>[x,y]
>>plot（x,y,'bo-'）
结果及其图像
ans=
00.5000
0.15000.7402
0.30001.0121
0.45001.3136
0.60001.6422
0.75001.9950
0.90002.3687
1.05002.7594
1.20003.1625
1.35003.5728
1.50003.9842
1.65004.3896
1.80004.7809
1.95005.1484
四．步长0.2
在MATLAB输入以下程序：
>>clearall
>>fun=inline（'y-x^2+1'）;
>>[x,y]=gjEuler（fun,[0,2],0.5,0.2）;
>>[x,y]
>>plot（x,y,'yd-'）
结果及其图像
ans=
00.5000
0.20000.8260
0.40001.2069
0.60001.6372
0.80002.1102
1.00002.6177
1.20003.1496
1.40003.6937
1.60004.2351
1.80004.7556
2.00005.2331
五．不同步长之比图像
六．由以上数值分析结果绘制表格
精确解
步长为0.05
步长为0.1
步长为0.15
步长为0.2
xi
yi
yi
误差
yi
误差
yi
误差
yi
误差
0.0000
0.5000
0.5000
0
0.5000
0.0000
0.5000
0.0000
0.5000
0.0000
0.0500
0.5769
0.5768
-0.0001
0.1000
0.6574
0.6573
-0.0001
0.6570
-0.0004
0.1500
0.7416
0.7414
-0.0002
0.7402
-0.0014
0.2000
0.8293
0.8291
-0.0002
0.8284
-0.0009
0.8260
-0.0033
0.2500
0.9205
0.9202
-0.0003
0.3000
1.0151
1.0147
-0.0004
1.0137
-0.0014
1.0121
-0.0030
0.3500
1.1130
1.1126
-0.0004
0.4000
1.2141
1.2136
-0.0005
1.2122
-0.0019
1.2069
-0.0072
0.4500
1.3183
1.3178
-0.0005
1.3136
-0.0047
0.5000
1.4256
1.4250
-0.0006
1.4232
-0.0024
0.5500
1.5359
1.5352
-0.0007
0.6000
1.6489
1.6482
-0.0007
1.6459
-0.0030
1.6422
-0.0067
1.6372
-0.0117
0.6500
1.7647
1.7639
-0.0008
0.7000
1.8831
1.8822
-0.0009
1.8794
-0.0037
0.7500
2.0040
2.0030
-0.001
1.9950
-0.0090
0.8000
2.1272
2.1261
-0.0011
2.1228
-0.0044
2.1102
-0.0170
0.8500
2.2527
2.2514
-0.0013
0.9000
2.3802
2.3789
-0.0013
2.3750
-0.0052
2.3687
-0.0115
0.9500
2.5096
2.5082
-0.0014
1.0000
2.6409
2.6393
-0.0016
2.6348
-0.0061
2.6177
-0.0232
1.0500
2.7737
2.7720
-0.0017
2.7594
-0.0143
1.1000
2.9079
2.9061
-0.0018
2.9010
-0.0069
1.1500
3.0434
3.0415
升级会员 