常微分方程的数值解法（欧拉法、改进欧拉法、泰勒方法和龙格-库塔法）.doc
《常微分方程的数值解法（欧拉法、改进欧拉法、泰勒方法和龙格-库塔法）.doc》由会员分享，可在线阅读，更多相关《常微分方程的数值解法（欧拉法、改进欧拉法、泰勒方法和龙格-库塔法）.doc（11页珍藏版）》请在冰点文库上搜索。
［例1］用欧拉方法与改进的欧拉方法求初值问题
在区间[0,1]上取的数值解。
［解］ 欧拉方法的计算公式为
使用excel表格进行运算，相应结果如下
例一:
欧拉法
n
x
y
精确解
0
0
1
1
1
0.1
1
1.003322
2
0.2
1.006667
1.013159
3
0.3
1.019824
1.029142
4
0.4
1.039054
1.050718
5
0.5
1.063754
1.077217
6
0.6
1.093211
1.107932
7
0.7
1.126681
1.142165
8
0.8
1.163443
1.179274
9
0.9
1.202845
1.218689
10
1
1.244314
1.259921
现用matlab编程，程序如下
x0=0;
y0=1;
x
（1）=0.1;
y
（1）=y0+0.1*2*x0/（3*y0^2）;
forn=1:
9
x（n+1）=0.1*（n+1）;
y（n+1）=y（n）+0.1*2*x（n）/（3*y（n）^2）;
end;
x
y
结果为
x=
Columns1through8
0.10000.20000.30000.40000.50000.60000.70000.8000
Columns9through10
0.90001.0000
y=
Columns1through8
1.00001.00671.01981.03911.06381.09321.12671.1634
Columns9through10
1.20281.2443
改进的欧拉方法其计算公式为
本题的精确解为，可用来检验数值解的精确度，列出计算结果。
使用excel表格进行运算，相应如下
例一:
改进的欧拉法
n
x
预测y
校正y
精确解
0
0
1
1
1
1
0.1
1
1.003333
1.003322
2
0.2
1.009956
1.01318
1.013159
3
0.3
1.026169
1.029171
1.029142
4
0.4
1.048054
1.050751
1.050718
5
0.5
1.074904
1.077252
1.077217
6
0.6
1.105976
1.107965
1.107932
7
0.7
1.140549
1.142194
1.142165
8
0.8
1.177965
1.179297
1.179274
9
0.9
1.217646
1.218706
1.218689
10
1
1.259103
1.25993
1.259921
现用matlab编程，程序如下
x0=0;
y0=1;
x
（1）=0.1;
ya
（1）=y0+0.1*2*x0/（3*y0^2）;
y
（1）=y0+0.05*（2*x0/（3*y0^2）+2*x0/（3*ya^2））;
forn=1:
9
x（n+1）=0.1*（n+1）;
ya（n+1）=ya（n）+0.1*2*x（n）/（3*ya（n）^2）;
y（n+1）=y（n）+0.05*（2*x（n）/（3*y（n）^2）+2*x（n+1）/（3*ya（n+1）^2））;
end;
x
y
结果为
x=
Columns1through8
0.10000.20000.30000.40000.50000.60000.70000.8000
Columns9through10
0.90001.0000
y=
Columns1through8
1.00001.00991.02611.04791.07481.10591.14071.1783
Columns9through10
1.21831.2600
［例2］用泰勒方法解
分别用二阶、四阶泰勒方法计算点＝0.1,0.2,…,1.0处的数值解，并与精确解进行比较。
解：
二阶泰勒方法
对于本题
故
使用excel表格进行运算，相应结果如下
n
x
y
精确解
0
0
1
1
1
0.1
1.003333
1.003322
2
0.2
1.013223
1.013159
3
0.3
1.029291
1.029142
4
0.4
1.050969
1.050718
5
0.5
1.077575
1.077217
6
0.6
1.108388
1.107932
7
0.7
1.142704
1.142165
8
0.8
1.179878
1.179274
9
0.9
1.21934
1.218689
10
1
1.260601
1.259921
现用matlab编程，程序如下
x0=0;
y0=1;
x
（1）=0.1;
y
（1）=y0+0.1/（3*y0^2）*（2*x0+0.1*（1-4*x0^2/（3*y0^3）））;
forn=1:
9
x（n+1）=0.1*（n+1）;
y（n+1）=y（n）+0.1/（3*y（n）^2）*（2*x（n）+0.1*（1-4*x（n）^2/（3*y（n）^3）））;
end;
x
y
结果为
x=
Columns1through9
0.10000.20000.30000.40000.50000.60000.70000.80000.9000
Column10
1.0000
y=
Columns1through9
1.00331.01321.02931.05101.07761.10841.14271.17991.2193
Column10
1.2606
四阶泰勒方法
使用excel表格进行运算，相应结果如下
例二：
四阶泰勒方法
n
x
y
y的一阶导数
y的二阶导数
y的三阶导数
y的四阶导数
精确解
0
0
1
0
0.666667
0
-1.33333
1
1
0.1
1.003328
0.066225
0.653509
-0.13402
-1.18306
1.003322
2
0.2
1.013191
0.129884
0.616121
-0.2711
-0.79041
1.013159
3
0.3
1.029211
0.188808
0.560087
-0.40991
-0.2947
1.029142
4
0.4
1.050823
0.241496
0.49274
-0.54379
0.159787
1.050718
5
0.5
1.077346
0.287189
0.421266
-0.66342
0.482809
1.077217
6
0.6
1.108063
0.325785
0.351405
-0.76055
0.651545
1.107932
7
0.7
1.142274
0.357656
0.286967
-0.83057
0.690167
1.142165
8
0.8
1.179339
0.383461
0.229962
-0.87296
0.641412
1.179274
9
0.9
1.218692
0.403983
0.181038
-0.8903
0.54637
1.218689
10
1
1.25985
0.420021
0.13996
-0.88694
0.435585
1.259921
现用matlab编程，程序如下
x0=0;
y0=1;
ya0=2*x0/（3*y0^2）;%%一阶导数
yb0=2/（3*y0^2）-8*x0^2/（9*y0^5）;%%二阶导数
yc0=-4*x0/（3*y0^5）-80*x0^3/（27*y0^8）;%%三阶导数
yd0=-4/（3*y0^5）+40*x0^2/（3*y0^8）-1280*x0^4/（81*y0^11）;%%四阶导数
x
（1）=0.1;
y
（1）=y0+0.1*ya0+0.01/2*yb0+0.001/6*yc0+0.0001/24*yd0;
ya
（1）=2*x
（1）/（3*y
（1）^2）;%%一阶导数
yb
（1）=2/（3*y
（1）^2）-8*x
（1）^2/（9*y
（1）^5）;%%二阶导数
yc
（1）=-4*x
（1）/（3*y
（1）^5）-80*x
（1）^3/（27*y
（1）^8）;%%三阶导数
yd
（1）=-4/（3*y
（1）^5）+40*x
（1）^2/（3*y
（1）^8）-1280*x
（1）^4/（81*y
（1）^11）;%%四阶导数
forn=1:
9
x（n+1）=0.1*（n+1）;
y（n+1）=y（n）+0.1*ya（n）+0.01/2*yb（n）+0.001/6*yc（n）+0.0001/24*yd（n）;
ya（n+1）=2*x（n+1）/（3*y（n+1）^2）;%%一阶导数
yb（n+1）=2/（3*y（n+1）^2）-8*x（n+1）^2/（9*y（n+1）^5）;%%二阶导数
yc（n+1）=-4*x（n+1）/（3*y（n+1）^5）-80*x（n+1）^3/（27*y（n+1）^8）;%%三阶导数
yd（n+1）=-4/（3*y（n+1）^5）+40*x（n+1）^2/（3*y（n+1）^8）-1280*x（n+1）^4/（81*y（n+1）^11）;%%四阶导数
end;
x
Y
结果为
x=
Columns1through8
0.10000.20000.30000.40000.50000.60000.70000.8000
Columns9through10
0.90001.0000
y=
Columns1through8
1.00331.01321.02921.05081.07731.10811.14231.1793
Columns9through10
1.21871.2598
［例3］用标准四阶R－K方法求
在区间[0,1]上，的数值解以及在区间[1,10]上，的数值解，并与精确解进行比较。
解：
对于本题
使用excel表格进行运算，相应结果如下
例三：
标准四阶R－K方法
n
x
y
k1
k2
k3
k4
精确解
0
0
1
0
0.003333
0.003322
0.006623
1
1
0.1
1.003322
0.006623
0.009869
0.009837
0.012989
1.003322
2
0.2
1.013159
0.012989
0.01603
0.015983
0.018883
1.013159
3
0.3
1.029143
0.018883
0.021632
0.021575
0.024154
1.029142
4
0.4
1.050718
0.024154
0.02656
0.0265
0.028726
1.050718
5
0.5
1.077217
0.028726
0.030772
0.030715
0.032586
1.077217
6
0.6
1.107932
0.032586
0.034286
0.034234
0.035772
1.107932
7
0.7
1.142165
0.035772
0.037155
0.037111
0.03835
1.142165
8
0.8
1.179274
0.03835
0.039454
0.039417
0.040398
1.179274
9
0.9
1.218689
0.040399
0.041264
0.041235
0.041997
1.218689
10
1
1.259921
0.041997
0.042663
0.042641
0.043222
1.259921
现用matlab编程，程序如下
x0=0;
y0=1;
k10=2*0.1*x0/（3*y0^2）;
k20=2*0.1*（x0+0.05）/（3*（y0+k10/2）^2）;
k30=2*0.1*（x0+0.05）/（3*（y0+k20/2）^2）;
k40=2*0.1*（x0+0.1）/（3*（y0+k30）^2）;
x
（1）=0.1;
y
（1）=y0+（k10+2*k20+2*k30+k40）/6;
k1
（1）=2*0.1*x
（1）/（3*y
（1）^2）;
k2
（1）=2*0.1*（x
（1）+0.05）/（3*（y
（1）+k1
（1）/2）^2）;
k3
（1）=2*0.1*（x
（1）+0.05）/（3*（y
（1）+k2
（1）/2）^2）;
k4
（1）=2*0.1*（x
（1）+0.1）/（3*（y
（1）+k3
（1））^2）;
forn=1:
9
x（n+1）=0.1*（n+1）;
y（n+1）=y（n）+（k1（n）+2*k2（n）+2*k3（n）+k4（n））/6;
k1（n+1）=2*0.1*x（n+1）/（3*y（n+1）^2）;
k2（n+1）=2*0.1*（x（n+1）+0.05）/（3*（y（n+1）+k1（n+1）/2）^2）;
k3（n+1）=2*0.1*（x（n+1）+0.05）/（3*（y（n+1）+k2（n+1）/2）^2）;
k4（n+1）=2*0.1*（x（n+1）+0.1）/（3*（y（n+1）+k3（n+1））^2）;
end;
x
y
结果为
x=
Columns1through8
0.10000.20000.30000.40000.50000.60000.70000.8000
Columns9through10
0.90001.0000
y=
Columns1through8
1.00331.01321.02911.05071.07721.10791.14221.1793
Columns9through10
1.21871.2599
您好，欢迎您阅读我的文章，本WORD文档可编辑修改，也可以直接打印。
阅读过后，希望您提出保贵的意见或建议。
阅读和学习是一种非常好的习惯，坚持下去，让我们共同进步。
升级会员 