matlab第二章__常微分方程的數(shù)值解法.ppt

1、微分方程的數(shù)值解法,四階龍格庫塔法 （The Fourth-Order RungeKutta Method）,常微分方程（Ordinary differential equations, ODE),初值問題-給出初始值 邊值問題-給出邊界條件,與初值常微分方程解算有關(guān)的指令 ode23 ode45 ode113 ode23t ode15s ode23s ode23tb,一.解ODE的基本機(jī)理:,2. 把高階方程轉(zhuǎn)換成一階微分方程組,1. 列出微分方程,初始條件,令,(2.1),(2.2),(2.3),例：著名的Van der Pol方程,令,降為一階,初始條件,3. 根據(jù)式(2.2)編寫計(jì)算導(dǎo)

2、數(shù)的M函數(shù)文件-ODE文件,把t,Y作為輸入宗量，把 作為輸出宗量,%M function file name: dYdt.m function Yd = f (t, Y) Yd = f (t,Y) 的展開式,例Van der Pol方程,%M function file name: dYdt.m function Yd = f (t, Y) Yd=zeros(size(Y);,4. 使編寫好的ODE函數(shù)文件和初值 供微分方程解算指令（solver）調(diào)用,Solver解算指令的使用格式,輸出宗量形式,說明： t0：初始時刻；tN：終點(diǎn)時刻Y0：初值； tol：計(jì)算精度,例題1：著名的Van d

3、er Pol方程,% 主程序 （程序名：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),解法1：采用ODE命令,Van der Pol方程,% 子程序 （程序名： dYdt.m ） function Ydot = dYdt (t, Y),Ydot=Y(2);-Y(2)*(Y(1)2-1)-Y(1);,或?qū)憺?fun

4、ction Ydot = dYdt (t, Y) Ydot=zeros(size(Y); Ydot(1)=Y(2); Ydot(2)=-Y(2)*(Y(1).2-1)-Y(1);,各種solver 解算指令的特點(diǎn),二. 四 階 Runge-Kutta 法,對 I=a,b作分割,步長,單步法-Runge-Kutta 方法 多步法-Admas方法,計(jì)算 的近似值 時只用到 ，是自開始方法,Runge-Kutta法是常微分方程的一種經(jīng)典解法 MATLAB 對應(yīng)命令：ode45,四階Runge-Kutta公式,四 階 Runge-Kutta 法計(jì)算流程圖,開始,Plot,初始條件： ； 積分步長： 迭

5、代次數(shù)：,輸出結(jié)果,子程序計(jì)算,End,三. Runge-Kutta 法解Van der Pol 方程的Matlab 程序結(jié)構(gòu)主程序：RK_vanderpol.m 子程序：RK_sub.m（函數(shù)文件）,解法2：采用Runge_Kutta法編程計(jì)算,主程序：RK_vanderpol.m t0=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); K2 = RK_sub(t0 + h/2, y0 + h*K1/2);

6、 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; end subplot (121), plot (t, yy1, t, yy2); grid subplot (122), plot (yy2, yy1); grid,子程序：RK_sub.m function ydot = vdpol (t, y) ydot=zeros(size

7、(y); ydot(1) = y(2); ydot(2) = -y(2)*(y(1)2-1)-y(1); 或?qū)憺椋?ydot = y(1) ;-y(2)*(y(1)2-1)-y(1);,四. Matlab對應(yīng)命令：ode23,ode45,調(diào)用格式： t, y=ode23 (函數(shù)文件名, t0, tN, y0, tol) t, y=ode45 (函數(shù)文件名, t0, tN, y0, tol),默認(rèn)精度： ode231e-3 ode451e-6,說明： t0：初始時刻；tN：終點(diǎn)時刻 y0：初值； tol：計(jì)算精度,3月15日作業(yè): 1.Van der Pol 方程的兩種解法:1)采用ode45命

8、令 2)Runge-Kutta方法 2.Duffing 方程的求解（Runge-Kutta方法，計(jì)算步長h=0.005，計(jì)算時間t0=0.0,tN=100） 要求：寫出程序體，打印所繪圖形，圖形標(biāo)題用個人的名字。,Duffing 方程,五. 動力學(xué)系統(tǒng)的求解,1. 動力學(xué)方程,2. 二階方程轉(zhuǎn)成一階方程,(1),令：,(2),其中：,即：,(2),3. Matlab 程序(主程序：ZCX),t0;Y0;h;N;P0,w; %輸入初始值、步長、迭代次數(shù)、初始激勵力； for i = 1 : N t1 = t0 + h P=P0*sin(w*t0);0.0;0.0 %輸入t0時刻的外部激勵力 K1

9、 = 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= %輸入 (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ù)或圖形,Matlab 程序(子程序：ZCX_sub.m),function ydot =

10、f (t, Y,P) M=, K=, C= %輸入結(jié)構(gòu)參數(shù) P1=zeros(3,1);inv（M）*P； A=zeros(0,0), eye(n,n) ; -M-1K, -M-1C ydot =AY+P1,例題2：三自由度質(zhì)量彈簧系統(tǒng),矩陣表示,其中：,動力學(xué)方程：,解析解：,已知參數(shù)：m1=m2=m3=1, k1=2, k2=2, K3=1, K4=2, P0=1, 要求：采用四階龍格庫塔法編程計(jì)算三個質(zhì)量的響應(yīng)時程.計(jì)算時間 0 50,例如：,4階龍格庫塔法的結(jié)果,ode45 的結(jié)果,第一個質(zhì)量的位移響應(yīng)時程,結(jié)果完全一致,MATLAB程序 （1）4階RK方法： （2）采用ode45：

11、m_chap2_ex2_1.m，m_chap2_ex2_1_sub.m,例題3： 蹦極跳系統(tǒng)的動態(tài)仿真,蹦極者系著一根彈性繩從高處的橋梁（或山崖等）向下跳。在下落的過程中，蹦極者幾乎處于失重狀態(tài)。按照牛頓運(yùn)動規(guī)律，自由下落的物體由下式確定：,其中，m 為人體的質(zhì)量，g 為重力加速度，x 為物體的位置，第二項(xiàng)和第三項(xiàng)表示空氣的阻力。其中位置 x 的基準(zhǔn)為橋梁的基準(zhǔn)面（即選擇橋梁作為位置的起點(diǎn) x0），低于橋梁的位置為正值，高于橋梁的位置為負(fù)值。如果人體系在一個彈性常數(shù)為 k 的彈性繩索上，定義繩索下端的初始位置為 0，則其對落體位置的影響為：,空氣的阻力,整個蹦極系統(tǒng)的數(shù)學(xué)模型為：,設(shè)橋梁距離地

12、面為 50 m，即 h2=50，蹦極者的起始位置為繩索的長度 30 m，即 h1=30，蹦極者起始速度為 0，其余的參數(shù)分別為 k20， a2a11；m70 kg，g10 m/s2。,初始條件：,已知參數(shù)：,初始條件變?yōu)椋?y0=-30; 0; % 初始位移和初始速度 t,y=ode45(bengji_sub, 0:0.01:100, y0); x1=50. - y(:,1); % x1代表蹦極者與地面之間的距離 plot(t,x1); grid plot(t,y(:,1); grid % y(:,1)代表位移,主程序 (程序名：bengji.m),Matlab程序,function ydot=f(t,y) m=70; k=20; a1=1; a2=1; g=10; x=y(1); % x代表蹦極者的位移 x_dot=y(2); % x_dot 代表 x 的速度 if x0 ydot=0,1;-k/m,-a1/m-(a2/m)*abs(x_dot)*y+0;g; else ydot=0,1;0,-a1/m-(a2/m)*abs(x_dot)*y+0;g; end,子程序 （程序名：b