常微分方程地数值解法欧拉法改进欧拉法泰勒方法和龙格库塔法Word格式.docx
《常微分方程地数值解法欧拉法改进欧拉法泰勒方法和龙格库塔法Word格式.docx》由会员分享,可在线阅读,更多相关《常微分方程地数值解法欧拉法改进欧拉法泰勒方法和龙格库塔法Word格式.docx(11页珍藏版)》请在冰豆网上搜索。
![常微分方程地数值解法欧拉法改进欧拉法泰勒方法和龙格库塔法Word格式.docx](https://file1.bdocx.com/fileroot1/2023-1/26/a14a4e25-59e3-470f-9925-cb010f8329ee/a14a4e25-59e3-470f-9925-cb010f8329ee1.gif)
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.244314
1.259921
现用matlab编程,程序如下
x0=0;
y0=1;
x
(1)=0.1;
y
(1)=y0+0.1*2*x0/(3*y0^2);
forn=1:
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=
1.00001.00671.01981.03911.06381.09321.12671.1634
1.20281.2443
改进的欧拉方法其计算公式为
本题的精确解为
,可用来检验数值解的精确度,列出计算结果。
使用excel表格进行运算,相应如下
改进的欧拉法
预测y
校正y
1.003333
1.009956
1.01318
1.026169
1.029171
1.048054
1.050751
1.074904
1.077252
1.105976
1.107965
1.140549
1.142194
1.177965
1.179297
1.217646
1.218706
1.259103
1.25993
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));
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));
1.00001.00991.02611.04791.07481.10591.14071.1783
1.21831.2600
[例2]用泰勒方法解
分别用二阶、四阶泰勒方法计算点
=0.1,0.2,…,1.0处的数值解,并与精确解进行比较。
解:
二阶泰勒方法
对于本题
故
1.013223
1.029291
1.050969
1.077575
1.108388
1.142704
1.179878
1.21934
1.260601
y
(1)=y0+0.1/(3*y0^2)*(2*x0+0.1*(1-4*x0^2/(3*y0^3)));
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)));
Columns1through9
0.10000.20000.30000.40000.50000.60000.70000.80000.9000
Column10
1.0000
1.00331.01321.02931.05101.07761.10841.14271.17991.2193
1.2606
四阶泰勒方法
例二:
y的一阶导数
y的二阶导数
y的三阶导数
y的四阶导数
0.666667
-1.33333
1.003328
0.066225
0.653509
-0.13402
-1.18306
1.013191
0.129884
0.616121
-0.2711
-0.79041
1.029211
0.188808
0.560087
-0.40991
-0.2947
1.050823
0.241496
0.49274
-0.54379
0.159787
1.077346
0.287189
0.421266
-0.66342
0.482809
1.108063
0.325785
0.351405
-0.76055
0.651545
1.142274
0.357656
0.286967
-0.83057
0.690167
1.179339
0.383461
0.229962
-0.87296
0.641412
1.218692
0.403983
0.181038
-0.8903
0.54637
1.25985
0.420021
0.13996
-0.88694
0.435585
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);
%%四阶导数
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);
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);
Y
1.00331.01321.02921.05081.07731.10811.14231.1793
1.21871.2598
[例3]用标准四阶R-K方法求
在区间[0,1]上,
的数值解以及在区间[1,10]上,
的数值解,并与精确解进行比较。
例三:
标准四阶R-K方法
k1
k2
k3
k4
0.003333
0.003322
0.006623
0.009869
0.009837
0.012989
0.01603
0.015983
0.018883
1.029143
0.021632
0.021575
0.024154
0.02656
0.0265
0.028726
0.030772
0.030715
0.032586
0.034286
0.034234
0.035772
0.037155
0.037111
0.03835
0.039454
0.039417
0.040398
0.040399
0.041264
0.041235
0.041997
0.042663
0.042641
0.043222
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);
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);
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);
1.00331.01321.02911.05071.07721.10791.14221.1793
1.21871.2599