微分方程的數(shù)值解法matlab(四階龍格—庫塔法)匯編

1、常微分方程常微分方程（Ordinary differential equations, ODE)n初值問題初值問題-給出初始值給出初始值n邊值問題邊值問題-給出邊界條件給出邊界條件與初值常微分方程解算有關(guān)的指令與初值常微分方程解算有關(guān)的指令ode23 ode45 ode113 ode23tode15s ode23sode23tb一一.解解ODE的基本機(jī)理的基本機(jī)理:2. 把高階方程轉(zhuǎn)換成一階微分方程組把高階方程轉(zhuǎn)換成一階微分方程組1. 列出微分方程列出微分方程初始條件初始條件令令)0()0(210yyYyyyy21,),(),(),(2121YtfYtfYtfyyY(2.1) (2.2)(2.

2、3),(tyyfy 00)0(,)0(yyyy例：著名的例：著名的Van der Pol方程方程0)1(2yyyy 令令 yyyy21,降為一階降為一階1221221) 1(yyyyyyY21yyY初始條件初始條件2010210)0()0(yyyyY3. 根據(jù)式根據(jù)式(2.2)編寫計算導(dǎo)數(shù)的編寫計算導(dǎo)數(shù)的M函數(shù)文件函數(shù)文件-ODE文件文件把把t,Y作為輸入宗量，把作為輸入宗量，把 作為輸出宗量作為輸出宗量Y %M function file name: dYdt.m function Yd = f (t, Y) Yd = f (t,Y) 的展開式的展開式例例Van der Pol方程方程 %M

3、 function file name: dYdt.m function Yd = f (t, Y) Yd=zeros(size(Y););1 () 2 (*) 12).1 () 2 ();2 () 1 (YYYYdYYd4. 使編寫好的使編寫好的ODE函數(shù)文件和初值函數(shù)文件和初值 供微分供微分方程解算指令（方程解算指令（solver）調(diào)用）調(diào)用0YSolver解算指令的使用格式解算指令的使用格式t, Y=solver (ODE函數(shù)文件名函數(shù)文件名, t0, tN, Y0, tol);ode45)()()()()()(222112110201tytytytytytyY輸出宗量形式輸出宗量形式)

4、() 1(:,1tyY)()2(:,2tyY說明：說明：t0：初始時刻；：初始時刻；tN：終點(diǎn)時刻：終點(diǎn)時刻Y0：初值；：初值； tol：計算精度：計算精度例題例題1：著名的著名的Van der Pol方程方程0)1(2yyyy % 主程序主程序 （程序名：程序名：VanderPol _ex1.m） t0 = 0; tN = 20; tol = 1e-6; Y0 = 0.25; 0.0; t, Y=ode45 (dYdt, t0, tN, Y0, tol); subplot (121), plot (t, Y) subplot (122), plot (Y( :, 1), Y( :, 2)解法

5、解法1：采用：采用ODE命令命令Van der Pol方程方程% 子程序子程序 （程序名：程序名： dYdt.m ） function Ydot = dYdt (t, Y)Ydot=Y(2);-Y(2)*(Y(1)2-1)-Y(1);或?qū)憺榛驅(qū)憺閒unction Ydot = dYdt (t, Y)Ydot=zeros(size(Y);Ydot(1)=Y(2);Ydot(2)=-Y(2)*(Y(1).2-1)-Y(1);解法指令解題類型特 點(diǎn)適合場合ode45非剛性非剛性采用采用4、5階階RungeKutta法法大多數(shù)場合的首選算法大多數(shù)場合的首選算法ode23非剛性非剛性采用采用Adams算

6、法算法較低精度（較低精度（103）場合）場合ode113非剛性非剛性多步法；采用多步法；采用Adams算法；高算法；高低精度均可（低精度均可（103106）ode45計算時間太長時計算時間太長時取代取代ode45ode23t適度剛適度剛性性采用梯形法則算法采用梯形法則算法適度剛性適度剛性ode15s剛性剛性多步法；采用多步法；采用2階階Rosenbrock算式，精度中等算式，精度中等當(dāng)當(dāng)ode45失敗時使用；失敗時使用；或存在質(zhì)量矩陣時或存在質(zhì)量矩陣時ode23s剛性剛性一步法；采用一步法；采用2階階Rosenbrock算式，低精度算式，低精度低精度時，比低精度時，比ode15s有有效；或存在

7、質(zhì)量矩陣時效；或存在質(zhì)量矩陣時ode23tb剛性剛性采用梯形法則反向數(shù)值微分采用梯形法則反向數(shù)值微分兩階段算法，低精度兩階段算法，低精度低精度時，比低精度時，比ode15s有有效；或存在質(zhì)量矩陣時效；或存在質(zhì)量矩陣時各種各種solver 解算指令的特點(diǎn)解算指令的特點(diǎn)二二. 四四 階階 Runge-Kutta 法法00)(),(ytybtaytfdtdybtttaN110對對 I=a,b作分割作分割1, 1 , 0,1Nihhhiii步長初值問題的數(shù)值初值問題的數(shù)值解法分為兩大類解法分為兩大類單步法-Runge-Kutta 方法多步法-Admas方法計算 的近似值 時只用到 ，是自開始方法 )(

8、1ntynnyt ,1nyuRunge-Kutta法是常微分方程的一種經(jīng)典解法常微分方程的一種經(jīng)典解法uMATLAB 對應(yīng)命令：對應(yīng)命令：ode45四階四階Runge-Kutta公式公式),()2,21()2,21(),()22(6342312143211hkyhtfkkhyhtfkkhyhtfkytfkkkkkhyynnnnnnnnnn四四 階階 Runge-Kutta 法計算流程圖法計算流程圖開始開始Next ifor i = 1 : N Plot)22(643211KKKKhyynn11,nnnnyytt初始條件：初始條件： ； 積分步長：積分步長： 迭代次數(shù)：迭代次數(shù)：0t0yhN輸出

9、結(jié)果httnn10ttn子程序計算iKEnd三三. Runge-Kutta 法解法解Van der Pol 方程的方程的Matlab 程序結(jié)構(gòu)程序結(jié)構(gòu)主程序：主程序：RK_vanderpol.m 子程序：子程序：RK_sub.m（函數(shù)文件）（函數(shù)文件） 解法解法2：采用：采用Runge_Kutta法編程計算法編程計算主程序：主程序：RK_vanderpol.mt0=0; tN=20; y0=0.25; 0; h=0.001;t = t0 : h : tN; N = length (t); j = 1;for i = 1 : N t1 = t0 + h; K1 = RK_sub(t0, y0);

10、 K2 = RK_sub(t0 + h/2, y0 + h*K1/2); K3 = RK_sub(t0 + h/2, y0 + h*K2/2); K4 = RK_sub(t0 + h, y0 + h*K3); y1 = y0 + (h/6)*(K1 + 2*K2 + 2*K3 + K4); yy1(j)=y1(1); yy2(j)=y1(2); t0=t1; y0=y1; j=j+1;endsubplot (121), plot (t, yy1, t, yy2); gridsubplot (122), plot (yy2, yy1); grid子程序：子程序：RK_sub.m function

11、 ydot = vdpol (t, y) ydot=zeros(size(y); ydot(1) = y(2); ydot(2) = -y(2)*(y(1)2-1)-y(1); 或?qū)憺椋夯驅(qū)憺椋?ydot = y(1) ;-y(2)*(y(1)2-1)-y(1);四四. Matlab對應(yīng)命令：對應(yīng)命令：ode23,ode45調(diào)用格式調(diào)用格式： t, y=ode23 (函數(shù)文件名函數(shù)文件名, t0, tN, y0, tol) t, y=ode45 (函數(shù)文件名函數(shù)文件名, t0, tN, y0, tol)默認(rèn)精度：默認(rèn)精度： ode231e-3 ode451e-6說明：說明：t0：初始時刻；：初

12、始時刻；tN：終點(diǎn)時刻：終點(diǎn)時刻y0：初值；：初值； tol：計算精度：計算精度3月月15日作業(yè)日作業(yè): 1.Van der Pol 方程的兩種解法方程的兩種解法:1)采用采用ode45命令命令 2)Runge-Kutta方法方法2.Duffing 方程的求解方程的求解（Runge-Kutta方法，計算步長方法，計算步長h=0.005，計算時間，計算時間t0=0.0,tN=100）要求：要求：寫出程序體，打印所繪圖形，圖形標(biāo)題用個寫出程序體，打印所繪圖形，圖形標(biāo)題用個人的名字。人的名字。)2.1cos(28.03.03tyyyy Duffing 方程方程0.0)0(,01.0)0(yy五五.

13、動力學(xué)系統(tǒng)的求解動力學(xué)系統(tǒng)的求解1. 動力學(xué)方程動力學(xué)方程2. 二階方程轉(zhuǎn)成一階方程二階方程轉(zhuǎn)成一階方程00)0(,)0()(),()()()(XXXXXPKXXCXM nRttttt(1) 0 () 0 () 0 (,)(,)()()(2XXYYXXYnRtttt令：令：)()()(tttPAYY(2)12112211)()(,nnnnnnnnRttRPM0PCMKMI0A其中：其中：)()()()()(1111tttttnnnnnPM0XXCMKMI0XX 即：即：00)0(,)0()(),()()()(XXXXXPKXXCXM nRttttt)()()(tttPAYY(2)3. Matl

14、ab 程序程序(主程序：主程序：ZCX) t0;Y0;h;N;P0,w; %輸入初始值、步長、迭代次數(shù)、初始激勵力；輸入初始值、步長、迭代次數(shù)、初始激勵力；for i = 1 : N t1 = t0 + h P=P0*sin(w*t0);0.0;0.0 %輸入輸入t0時刻的外部激勵力時刻的外部激勵力 K1 = ZCX_sub (t0, Y0, P ) P= %輸入輸入 (t0+h/2) 時刻的外部激勵力時刻的外部激勵力 K2 = ZCX_sub (t0 + h/2, Y0 + hK1/2, P ) K3 = ZCX_sub (t0 + h/2, Y0 + hK2/2, P ) P= %輸入輸入

15、 (t0+h) 時刻的外部激勵力時刻的外部激勵力 K4 = ZCX_sub (t0 + h, y0 + hK3, P) Y1 = y0 + (h/6) (K1 + 2K2 + 2K3 + K4) t1, Y1 (輸出輸出 t1, y1)next i輸出數(shù)據(jù)或圖形輸出數(shù)據(jù)或圖形)()()(tttPAYYMatlab 程序程序(子程序：子程序：ZCX_sub.m)function ydot = f (t, Y,P) M=, K=, C= %輸入結(jié)構(gòu)參數(shù)輸入結(jié)構(gòu)參數(shù) P1=zeros(3,1);inv（M）*P； A=zeros(0,0), eye(n,n) ; -M-1K, -M-1C ydot

16、=AY+P1)()()(tttPAYY例題例題2：三自由：三自由度質(zhì)量彈簧系統(tǒng)度質(zhì)量彈簧系統(tǒng)m1m2m3k1k2k3x1x2x3k4P0sin(wt)0)sin()(02121111tPxxkxkxm 0)()(32321222xxkxxkxm 0)(3432333xkxxkxm 00)sin(0000000003214333322221321321tPxxxkkkkkkkkkkxxxmmm 矩陣表示矩陣表示其中：其中：321xxxX321000000mmmM433332222100kkkkkkkkkkK00)sin()(0tPtP動力學(xué)方程：動力學(xué)方程：)sin(22221. 0)sin(3

17、0363. 2)sin(121088. 0)()()(000321twkPtwkPtwkPtxtxtx解析解：已知參數(shù)已知參數(shù)：m1=m2=m3=1, k1=2, k2=2, K3=1, K4=2, P0=1, 要求：要求：采用四階龍格庫塔法編程計算三個質(zhì)量的響應(yīng)時程.計算時間 0 50例如：例如：)()()(tttPKXXM 305101520253035404550-0.4-0.3-0.2-0.100.10.20.30.405101520253035404550-0.4-0.3-0.2-0.100.10.20.30.44階龍格庫塔法的結(jié)果階龍格庫塔法的結(jié)果ode45 的結(jié)果的結(jié)果第一個質(zhì)量

18、的位移響應(yīng)時程第一個質(zhì)量的位移響應(yīng)時程結(jié)果完全一致結(jié)果完全一致MATLAB程序程序（1）4階階RK方法：方法： （2）采用）采用ode45： m_chap2_ex2_1.m，m_chap2_ex2_1_sub.m 例題例題3： 蹦極跳系統(tǒng)的動態(tài)仿真蹦極跳系統(tǒng)的動態(tài)仿真蹦極者系著一根彈性繩從高處的橋梁（或山崖等）向下跳。在蹦極者系著一根彈性繩從高處的橋梁（或山崖等）向下跳。在下落的過程中，蹦極者幾乎處于失重狀態(tài)。按照牛頓運(yùn)動規(guī)律，下落的過程中，蹦極者幾乎處于失重狀態(tài)。按照牛頓運(yùn)動規(guī)律，自由下落的物體由下式確定自由下落的物體由下式確定：xxaxamgxm 210 , 00,)(xxkxxb其中，其

19、中，m 為人體的質(zhì)量，為人體的質(zhì)量，g 為重力加速度，為重力加速度，x 為物體的位置，第二項和第三項表示空氣的阻為物體的位置，第二項和第三項表示空氣的阻力。其中位置力。其中位置 x 的基準(zhǔn)為橋梁的基準(zhǔn)面（即選的基準(zhǔn)為橋梁的基準(zhǔn)面（即選擇橋梁作為位置的起點(diǎn)擇橋梁作為位置的起點(diǎn) x0），低于橋梁的），低于橋梁的位置為正值，高于橋梁的位置為負(fù)值。如果人位置為正值，高于橋梁的位置為負(fù)值。如果人體系在一個彈性常數(shù)為體系在一個彈性常數(shù)為 k 的彈性繩索上，定義的彈性繩索上，定義繩索下端的初始位置為繩索下端的初始位置為 0，則其對落體位置的，則其對落體位置的影響為：影響為：地面地面x橋梁基準(zhǔn)面橋梁基準(zhǔn)面0梯

20、梯子子h2h1空氣的阻力空氣的阻力xxaxaxbmgxm 21)(整個蹦極系統(tǒng)的數(shù)學(xué)模型為：整個蹦極系統(tǒng)的數(shù)學(xué)模型為：設(shè)橋梁距離地面為設(shè)橋梁距離地面為 50 m，即，即 h2=50，蹦極者的起始位置為繩索，蹦極者的起始位置為繩索的長度的長度 30 m，即，即 h1=30，蹦極者起始速度為，蹦極者起始速度為 0，其余的參數(shù)分，其余的參數(shù)分別為別為 k20， a2a11；m70 kg，g10 m/s2。地面地面x橋梁基準(zhǔn)面橋梁基準(zhǔn)面0梯梯子子h2h1初始條件：初始條件：0)0(;30)0(xx已知參數(shù)：已知參數(shù)：10,70, 1,20,50,302121gmaakhh0 , 00,)(xxkxxbxxY0,00100,0102121xgYxmamaYxgYxmamamkY令令： 0300)0()0(xxY初始條件變?yōu)椋撼跏紬l件變?yōu)椋簒xaxaxbmgxm 21)(0 , 00,)(xxkxxby0=-30; 0; % 初始位移和初始速度初始位移和初始速度t,y=ode45(bengji_sub, 0:0.01:100, y0);x1=50. - y(:,1); % x1代表蹦極者與地面之間的距離代表蹦極者與