常微分方程的数值解法及仿真.pdf-资料库

常微分方程的数值解法及仿真1 支强 zhistar@163.com 西安交通大学自动化系 2010 年 1 月 7 日 1 ©本文仅限用于学习交流，转载请注明出处。文中如有不正之处，欢迎指正！

西安交通大学科技报告 zhistar@163.com 目录一、欧拉（Euler）公式..................................................................2 二、龙格-库塔公式 .........................................................................2 1. 二阶龙格-库塔公式 .................................................................2 2. 四阶龙格-库塔公式 .................................................................2 三、一阶常微分方程组的数值解法 ..............................................2 四、仿真算例 ..................................................................................4 仿真 1 应用欧拉法 .......................................................................4 仿真 2 应用二阶龙格-库塔法......................................................5 仿真 3 应用四阶龙格-库塔法......................................................6 附录 Matlab程序 .................................................................................7 1. 欧拉法程序 ..............................................................................7 2. 二阶龙格-库塔法程序 .............................................................8 3. 四阶龙格-库塔法程序 .............................................................9 参考文献............................................................................................ 10 1

西安交通大学科技报告 zhistar@163.com 设常微分方程为 = dy dt y a ( ) ⎧ ⎪ ⎨ ⎪ ⎩ f t y ( , ) = y 0 一、欧拉（Euler）公式 Equation Chapter 1 Section 1 y n 1 + = y n + hf t y ( , n n ) (1.1) 二、龙格-库塔公式 1. 二阶龙格-库塔公式 Equation Section (Next) = f f y ⎧ ⎪ n 1 + ⎪ k =⎨ 1 ⎪ k = ⎪ ⎩ 2 2. 四阶龙格-库塔公式 y + n t y ( , n t ( + n k 1 + k ) 2 ( h 2 ) n h y , + hk 1 ) n (2.1) k 1 + 2 k 2 + 2 k 3 + k 4) y ⎧ ⎪ n 1 + ⎪ k =⎪ 1 ⎪⎪ k = ⎨ ⎪ ⎪ = k ⎪ ⎪ k = ⎪⎩ 2 3 4 = f y + n t y ( n , f t ( n + f f t ( n t ( n + + ( h 6 ) n h 2 h 2 h y , , , h 2 h 2 hk y n + y n + + n k 1 ) (2.2) k ) 2 ) 3 三、一阶常微分方程组的数值解法 Equation Section (Next) 将由个一阶组成的常微分方程初值问题： m 2

西安交通大学科技报告 zhistar@163.com dy 1 dt dy 2 dt dy m dt y a ( ) 1 y a ( ) 2 y a ( ) m ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ = f 1 t y y ( , , 1 2 , , y ) m = = = = = f 2 t y y ( , , 1 2 , , y ) m f m t y y ( , , 1 2 , , y ) m (3.1) η 1 η 2 η m ⎧ ⎪ ⎨ ⎪ ⎩ ⎡ ⎢ ⎢ ⎢ ⎢ ⎣ Y x ( ) = = ( , F t y dY dt Y a η ( ) = ) (3.2) 1 f x y y ( , , 1 2 f x y y , ( , 2 x y y ( , , 1 f 2 m 1 , , y m y m ) ) , y ) m 2 ， = η ⎤ ⎥ ⎥ ⎥ ⎥ ⎦ η ⎡ ⎤ 1 ⎢ ⎥ η ⎢ ⎥ 2 。 ⎢ ⎥ ⎢ ⎥ η ⎣ ⎦ m 写成向量形式：其中： Y x ( ) = ⎡ ⎢ ⎢ ⎢ ⎢ ⎣ y x ( ) 1 y x ( ) 2 ， x ( ) ⎤ ⎥ ⎥ ⎥ ⎥ ⎦ y m 下面以具有两个方程的方程组为例，给出欧拉方法、龙格-库塔公式。常微分方程组： ⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ = dy dt dz dt y a ( ) z a ( ) = f t y z ( , , ) g t y z ( , , ) = = y 0 z 0 1．欧拉（Euler）公式： y z ⎛ ⎜ ⎝ n 1 + n 1 + 2．预估-校正公式： ⎛ ⎜ ⎝ n 1 + y z n ⎛ ⎜ ⎝ 1 + + y z n 1 + n 1 + ⎞ ⎟ ⎠ = ⎛ ⎜ ⎝ y z n n ⎞ ⎟ ⎠ ⎞ ⎟ ⎠ ⎞ ⎟ ⎠ h 2 = = ⎛ ⎜ ⎝ ⎛ ⎜ ⎝ y z n n y z n n ⎞ ⎟ ⎠ ⎞ ⎟ ⎠ + h + h ⎛ ⎜ ⎝ ⎛ ⎜ ⎝ f t y z ( , , n g t y z , , ( n n n n n f t y z ( , , n g t y z , , ( n n n n n ) ) ⎞ ⎟ ⎠ ) ) ⎞ ⎟ ⎠ ⎡ ⎢ ⎣ ⎛ ⎜ ⎝ f t ( g t ( n n , , y z , n y z , n n n ) ) ⎞ ⎟ ⎠ + ⎛ ⎜ ⎝ f t ( g t ( n n , , y y , , z z n 1 + n 1 + n 1 + n 1 + 3．四阶龙格-库塔公式： (3.3) (3.4) (3.5) (3.6) 3 ) ) ⎞ ⎟ ⎠ ⎤ ⎥ ⎦

西安交通大学科技报告 zhistar@163.com K 1 + 2 K + 2 K 3 2 ] + 4 K [ h 6 k (1) 1 k (2) 1 ⎞ ⎟ ⎟ ⎠ + 2 ⎛ ⎜ ⎜ ⎝ k k (1) 2 (2) 2 ⎞ ⎟ ⎟ ⎠ + 2 Y k + = 1 Y n + ⎞ ⎟ ⎠ = ⎛ ⎜ ⎝ y z n n + ⎞ ⎟ ⎠ h 6 ⎡ ⎢ ⎢ ⎣ ⎛ ⎜ ⎜ ⎝ y z n 1 + n 1 + ⎛ ⎜ ⎜ ⎝ ) ) k (1) 3 k (2) 3 ⎞ ⎟ ⎟ ⎠ + ⎛ ⎜ ⎜ ⎝ k k (1) 4 (2) 4 ⎞ ⎟ ⎟ ⎠ ⎤ ⎥ ⎥ ⎦ ⎞ ⎟ ⎠ K 1 = ⎛ ⎜ ⎜ ⎝ k (1) 1 k (2) 1 ⎞ ⎟ ⎟ ⎠ = ⎛ ⎜ ⎝ f t ( g t ( n n , , y z , n y z , n n n K 2 = ⎛ ⎜ ⎜ ⎝ k k (1) 2 (2) 2 ⎞ ⎟ ⎟ ⎠ K 3 = ⎛ ⎜ ⎜ ⎝ k (1) 3 k (2) 3 ⎞ ⎟ ⎟ ⎠ ⎛ ⎜ = ⎜ ⎜ ⎜ ⎝ ⎛ ⎜ = ⎜ ⎜ ⎜ ⎝ f t ( n + g t ( n + f t ( n + g t ( n + h 2 h 2 h 2 h 2 , y n + , y n + , y n + , y n + h 2 h 2 h 2 h 2 k (1) 1 , z n + k (1) 1 , z n + k (1) 2 , z n + k (1) 2 , z n + h 2 h 2 h 2 h 2 k (2) 1 ) k (2) 1 ) k (2) 2 ) k (2) 2 ) ⎞ ⎟ ⎟ ⎟ ⎟ ⎠ ⎞ ⎟ ⎟ ⎟ ⎟ ⎠ K 4 = ⎛ ⎜ ⎜ ⎝ k k (1) 4 (2) 4 ⎞ ⎟ ⎟ ⎠ ⎛ = ⎜ ⎜ ⎝ t f ( g t ( n n + + h y , h y , n n + + hk hk (1) 3 (1) 3 , , z z n n + + hk hk (2) 3 (2) 3 ) ) ⎞ ⎟ ⎟ ⎠ ⎛ ⎜ ⎝ (3.7) (3.8) (3.9) (3.10) (3.11) (3.12) 四、仿真算例设有一过程 y xω ω= 2 − y ，其中， 20ω= ， (0) y = ， (0) 0 y = 0 。当输入 x为阶跃输入时，系统的理论输出为 Y theory = 1 (1 cos ω − tω ) 。本例出处见文献[1]第 5 页。 N ∑ n 1 = Y ( theory − Y simulation 2 ) N ,此式将作为评价数值积分算法性能的标准之定义均方误差 E = 一。仿真 1 应用欧拉法应用（1.1）式的仿真结果见图 4.1，E=1.2664e-007。 4

0.12 0.1 0.08 0.06 0.04 0.02 0 -0.02 0 西安交通大学科技报告 zhistar@163.com 欧拉法 Theory Practice 100 200 300 400 500 600 700 800 900 1000 图 4.1 向前欧拉公式仿真 2 应用二阶龙格-库塔法应用（2.1）式的仿真结果见图 4.2，E=4.3717e-011。可见，二阶龙格-库塔法较欧拉法精度提高了四个数量级。由于仿真采用的算例为二阶常微分方程，可分解为两个一阶常微分方程，令 1x y= ， 2x y= ，则有 x 2 2 ω ω x =⎧ 1 ⎨ x = ⎩ 2 编程时这两个常微分方程的安排次序对精度有很大影响。一步预测时会引入估计误差， (1) (2) x − y 先计算（1）式后计算（2）式，即先求后求，避免不准确的值进一步影响值；修 y y y y 正时先修正，然后利用修正后的修正。三阶、四阶龙格-库塔算法编程时也应该注意 y y y 这一点。 5

0.1 0.09 0.08 0.07 0.06 0.05 0.04 0.03 0.02 0.01 0 0 西安交通大学科技报告 zhistar@163.com 二阶龙格-库塔法 Theory Practice 100 200 300 400 500 600 700 800 900 1000 图 4.2 二阶龙格-库塔法仿真 3 应用四阶龙格-库塔法应用（2.2）式的仿真结果见图 4.3，E=2.7926e-019。可见，四阶龙格-库塔法较二阶龙格-库塔法精度上又有了非常大的提高。 6

0.1 0.09 0.08 0.07 0.06 0.05 0.04 0.03 0.02 0.01 0 0 西安交通大学科技报告 zhistar@163.com 四阶龙格-库塔法 Theory Simulation 100 200 300 400 500 600 700 800 900 1000 图 4.3 四阶龙格-库塔法附录 Matlab 程序 1. 欧拉法程序 %欧拉法求解常微分方程 %DDY=f(Y,x)=omega*x-omega^2*Y 二阶常微分方程，x 此处为阶跃输入，Y 为状态 %==================================================================== clear all clear clc omega=20; T=0; Y=0; YD=0; X=1; H=0.001; ErrorSum=0; n=1; while T<=1 YTheory(n)=(1-cos(omega*T))/omega; %理论输出值 7

资料库

常微分方程的数值解法及仿真.pdf

相关推荐

课程资源

热门标签

最新资料