加载中…
常微分方程初值问题的数值解法
(2015-06-15 22:37:30)
标签：
杂谈 |
-
课题名称
通过编写程序使用数值方法解常微分方程初值问题。
-
解决的问题
利用Euler方法和改进的Euler方法求解初值问题：
分别取步长0.1， 0.2， 0.4. 此初值问题的精确解为:
-
采用的数值方法
-
Euler方法
本实验采用了Euler方法和改进Euler方法。
对于一阶常微分方程初值问题：
Euler方法是解常微分方程初值问题最简单最古老的一种数值方法，其基本思路就是把上式中的导数项用差商逼近，从而将一个微分方程转化为一个代数方程，以便求解。
设在中取等距离结点h，在结点上，由上式有：
http://araleii.com/wp-content/uploads/2015/06/061515_1436_1.png
又由差商定义可得：
所以有：
http://araleii.com/wp-content/uploads/2015/06/061515_1436_2.png
用的近似值带入上式，则有计算的Euler公式：
http://araleii.com/wp-content/uploads/2015/06/061515_1436_3.png
通过对Euler法公式泰勒展开进行分析可以发现Euler方法为1阶方法。
-
改进的Euler方法
利用梯形公式对Euler方法进行改进，可得：
http://araleii.com/wp-content/uploads/2015/06/061515_1436_4.png
改进的Euler方法是用Euler方法先求一个预测值，再用这个预测值来计算，即：
http://araleii.com/wp-content/uploads/2015/06/061515_1436_5.png
用泰勒展开容易知道改进的Euler方法具有二阶精度。
-
-
算法程序
要实现的方法是Euler法和改进Euler法，源程序在实验报告附件中，这里展示部分核心代码。
首先是迭代式的实现函数，实现代码如下：
http://araleii.com/wp-content/uploads/2015/06/061515_1436_6.png
以及一个求解精确值的函数，用于比较和计算误差：
http://araleii.com/wp-content/uploads/2015/06/061515_1436_7.png
接下来是Euler法的实现过程，利用之前的式子：
http://araleii.com/wp-content/uploads/2015/06/061515_1436_8.png
然后是改进Euler法的实现过程，同样是利用之前的式子：
http://araleii.com/wp-content/uploads/2015/06/061515_1436_9.png
完成了以上工作之后，只需要编写一个简单的测试函数就可以调用了：
http://araleii.com/wp-content/uploads/2015/06/061515_1436_10.png
-
数值结果
对于初值问题：
实验测试结果如下：
http://araleii.com/wp-content/uploads/2015/06/061515_1436_11.png
分别去不同步长数值结果统计如下：
表格 1 h=0.1
x
精确解
Euler法解
改进Euler法解
Euler法误差
改进Euler法误差
0
3
3
3
0
0
0.1
2.991697
3
2.991667
0.008303
0.00003
0.2
2.967145
2.983333
2.967087
0.016189
0.000057
0.3
2.927397
2.950482
2.927322
0.023085
0.000076
0.4
2.874147
2.902639
2.874072
0.028493
0.000075
0.5
2.809627
2.841656
2.809583
0.032029
0.000044
0.6
2.736491
2.769955
2.736518
0.033464
0.000027
0.7
2.657655
2.690401
2.657804
0.032746
0.000148
0.8
2.576133
2.606147
2.576456
0.030014
0.000323
0.9
2.494853
2.520442
2.495402
0.025589
0.000549
1
2.416484
2.436434
2.417301
0.01995
0.000816
1.1
2.343285
2.356965
2.34439
0.01368
0.001105
1.2
2.27698
2.28438
2.278371
0.007399
0.00139
1.3
2.218693
2.220377
2.220337
0.001684
0.001645
1.4
2.168938
2.165922
2.170779
0.003016
0.001841
1.5
2.127674
2.121244
2.129634
0.00643
0.001961
1.6
2.094403
2.08591
2.096396
0.008493
0.001993
1.7
2.068304
2.058985
2.070244
0.009319
0.00194
1.8
2.04837
2.039217
2.050183
0.009153
0.001813
1.9
2.033534
2.025235
2.035165
0.008299
0.001632
2
2.022765
2.015705
2.024182
0.00706
0.001417
表格 2 h=0.2
x
精确解
Euler法解
改进Euler法解
Euler法误差
改进Euler法误差
0
3
3
3
0
0
0.2
2.96715
3
2.966667
0.032855
0.000478
0.4
2.87415
2.933333
2.873358
0.059187
0.000789
0.6
2.73649
2.807758
2.735935
0.071267
0.000556
0.8
2.57613
2.641782
2.576739
0.065648
0.000606
1
2.41648
2.461357
2.419284
0.044873
0.0028
1.2
2.27698
2.29411
2.282579
0.01713
0.005599
1.4
2.16894
2.161986
2.177031
0.006952
0.008093
1.6
2.0944
2.074672
2.103778
0.019731
0.009375
1.8
2.04837
2.027742
2.057483
0.020628
0.009113
2
2.02277
2.007904
2.030436
0.014861
0.007671
表格 3 h=0.4
x
精确解
Euler法解
改进Euler法解
Euler法误差
改进Euler法误差
0
3
3
3
0
0
0.4
2.874147
3
2.866667
0.125853
0.00748
0.8
2.576133
2.733333
2.57119
0.1572
0.004943
1.2
2.27698
2.326959
2.296998
0.049979
0.020018
1.6
2.094403
2.03513
2.1444
0.059273
0.049997
2
2.022765
1.990552
2.082701
0.032213
0.059936
表格 4 h=0.04
x
精确解
Euler法解
改进Euler法解
Euler法误差
改进Euler法误差
0
3
3
3
0
0
0.04
2.99867
3
2.998667
0.001333
0.000001
0.08
2.99468
2.997333
2.994677
0.002654
0.000002
0.12
2.98806
2.992012
2.98806
0.00395
0.000002
0.16
2.97886
2.984068
2.97886
0.005205
0.000003
0.2
2.96715
2.973549
2.967141
0.006404
0.000003
0.24
2.95299
2.960522
2.952984
0.007534
0.000004
0.28
2.93649
2.945071
2.936486
0.008582
0.000004
0.32
2.91776
2.927299
2.91776
0.009535
0.000003
0.36
2.89694
2.90732
2.896933
0.010384
0.000002
0.4
2.87415
2.885266
2.874146
0.01112
0.000001
0.44
2.84955
2.861284
2.849552
0.011733
0.000001
0.48
2.82331
2.835529
2.823314
0.012219
0.000004
0.52
2.7956
2.808172
2.795606
0.012574
0.000008
0.56
2.7666
2.77939
2.766609
0.012793
0.000013
0.6
2.73649
2.749369
2.736509
0.012878
0.000019
0.64
2.70547
2.718301
2.705498
0.012829
0.000026
0.68
2.67373
2.686383
2.673768
0.01265
0.000034
0.72
2.64147
2.653814
2.641513
0.012344
0.000043
0.76
2.60887
2.620794
2.608927
0.01192
0.000054
0.8
2.57613
2.58752
2.576199
0.011386
0.000066
0.84
2.54344
2.554187
2.543513
0.010752
0.000079
0.88
2.51096
2.520986
2.511049
0.01003
0.000092
0.92
2.47887
2.488098
2.478974
0.009231
0.000107
0.96
2.44733
2.455698
2.44745
0.008371
0.000123
1
2.41648
2.423948
2.416623
0.007463
0.000139
1.04
2.38648
2.392998
2.38663
0.006523
0.000155
1.08
2.35742
2.362985
2.357592
0.005565
0.000172
1.12
2.32943
2.334032
2.329616
0.004604
0.000188
1.16
2.30259
2.306245
2.302793
0.003655
0.000204
1.2
2.27698
2.279712
2.2772
0.002732
0.000219
1.24
2.25266
2.254507
2.252895
0.001846
0.000234
1.28
2.22968
2.230685
2.229923
0.001009
0.000247
1.32
2.20805
2.208284
2.208312
0.000231
0.000259
1.36
2.18781
2.187327
2.188077
0.000481
0.000269
1.4
2.16894
2.167818
2.169216
0.00112
0.000278
1.44
2.15143
2.14975
2.151717
0.001682
0.000284
1.48
2.13527
2.1331
2.135555
0.002166
0.000289
1.52
2.1204
2.117832
2.120694
0.00257
0.000292
1.56
2.1068
2.103903
2.107091
0.002896
0.000292
1.6
2.0944
2.091256
2.094693
0.003147
0.000291
1.64
2.08316
2.07983
2.083443
0.003327
0.000287
1.68
2.073
2.069557
2.073278
0.00344
0.000282
1.72
2.06386
2.060366
2.064133
0.003492
0.000275
1.76
2.05567
2.052181
2.055938
0.003491
0.000266
1.8
2.04837
2.044927
2.048626
0.003443
0.000256
1.84
2.04188
2.038529
2.042128
0.003354
0.000245
1.88
2.03614
2.032911
2.036377
0.003233
0.000233
1.92
2.03109
2.028001
2.031307
0.003085
0.000221
1.96
2.02665
2.02373
2.026854
0.002917
0.000208
2
2.02277
2.020031
2.022959
0.002734
0.000194
-
讨论和分析
将之前的数值结果进行统计和整理如下：
http://araleii.com/wp-content/uploads/2015/06/061515_1436_12.png
图 1 h=0.1解的比较
http://araleii.com/wp-content/uploads/2015/06/061515_1436_13.png
图 2 h=0.1误差比较
http://araleii.com/wp-content/uploads/2015/06/061515_1436_14.png
图 3 h=0.2解的比较
http://araleii.com/wp-content/uploads/2015/06/061515_1436_15.png
图 4 h=0.2误差比较
http://araleii.com/wp-content/uploads/2015/06/061515_1436_16.png
图 5 h=0.4解的比较
http://araleii.com/wp-content/uploads/2015/06/061515_1436_17.png
图 6 h=0.4误差比较
http://araleii.com/wp-content/uploads/2015/06/061515_1436_18.png
图 7 h=0.04解的比较
http://araleii.com/wp-content/uploads/2015/06/061515_1436_19.png
图 8 h=0.04误差比较
从统计结果中可以看出，两种方法的都能很好的逼近精确解，改进Euler方法比Euler方法确实误差小很多，之所以Euler方法在大约1.3之后误差又再次增加是因为误差又正变为了负，取绝对误差就体现出了误差再次增大的效果。下面再对课本上的例子进行测试：
统计结果如下：
表格 5 测试2，h=0.1
x
精确解
Euler法解
改进Euler法解
Euler法误差
改进Euler法误差
0
1
1
1
0
0
0.1
1.095445
1.1
1.09
0.004555
0.005445
0.2
1.183216
1.191818
1.173191
0.008602
0.010025
0.3
1.264911
1.277438
1.251067
0.012527
0.013844
0.4
1.341641
1.358213
1.324566
0.016572
0.017075
0.5
1.414214
1.435133
1.394305
0.020919
0.019909
0.6
1.48324
1.508966
1.460691
0.025727
0.022549
0.7
1.549193
1.580338
1.523984
0.031145
0.02521
0.8
1.612452
1.649783
1.584335
0.037332
0.028117
0.9
1.67332
1.717779
1.641804
0.044459
0.031516
1
1.732051
1.784771
1.696376
0.05272
0.035675
作图分析可见：
http://araleii.com/wp-content/uploads/2015/06/061515_1436_20.png
图 9 方程2，h=0.1解的比较
http://araleii.com/wp-content/uploads/2015/06/061515_1436_21.png
图 10 方程2，h=0.1误差比较
可以看到，Euler法和改进Euler法能够较好的得到逼近精确解的数值结果。随着迭代次数的增加，改进Euler方法确实能够更好的逼近精确解。
-
总结
本次实验通过对Euler法的改进Euler法的实现，实现了通过数值方法对常微分方程的解析解的逼近。通过实验可以看出，拥有2阶精度的改进Euler方法确实比朴素的Euler方法效果更好的，误差更小。实现的过程中，最开始时出现了小错误是因为我把x增加的语句提前了，导致每次用的x提前了一个步长，使得结果不够理想。写程序的时候还是要仔细。然后通过计算机求解近似解难免会存在浮点误差，不过这往往是随着步长的调整一般能够得到解决的。
喜欢
0
赠金笔