微分方程的數值解法matlab(四階龍格—庫塔法)PPT課件

微分方程的數值解法 四階龍格 庫塔法 TheFourth OrderRunge KuttaMethod 常微分方程 Ordinarydifferentialequations ODE 初值問題 給出初始值邊值問題 給出邊界條件 與初值常微分方程解算有關的指令ode23ode45ode113ode23tode15sode23sode23tb 一 解ODE的基本機理 2 把高階方程轉換成一階微分方程組 1 列出微分方程 初始條件 令 2 1 2 2 2 3 例 著名的VanderPol方程 令 降為一階 初始條件 3 根據式 2 2 編寫計算導數的M函數文件 ODE文件 把t Y作為輸入宗量 把作為輸出宗量 Mfunction dYdt mfunctionYd f t Y Yd f t Y 的展開式 例VanderPol方程 Mfunction dYdt mfunctionYd f t Y Yd zeros size Y 4 使編寫好的ODE函數文件和初值供微分方程解算指令 solver 調用 Solver解算指令的使用格式 輸出宗量形式 說明 t0 初始時刻 tN 終點時刻Y0 初值 tol 計算精度 例題1 著名的VanderPol方程 主程序 程序名 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命令 VanderPol方程 子程序 程序名 dYdt m functionYdot dYdt t Y Ydot Y 2 Y 2 Y 1 2 1 Y 1 或寫為 functionYdot dYdt t Y Ydot zeros size Y Ydot 1 Y 2 Ydot 2 Y 2 Y 1 2 1 Y 1 各種solver解算指令的特點 二 四階Runge Kutta法 對I a b 作分割 步長 單步法 Runge Kutta方法多步法 Admas方法 計算的近似值時只用到 是自開始方法 Runge Kutta法是常微分方程的一種經典解法MATLAB對應命令 ode45 四階Runge Kutta公式 四階Runge Kutta法計算流程圖 開始 Plot 初始條件 積分步長 迭代次數 輸出結果 子程序計算 End 三 Runge Kutta法解VanderPol方程的Matlab程序結構主程序 RK vanderpol m子程序 RK sub m 函數文件 解法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 fori 1 Nt1 t0 h K1 RK sub t0 y0 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 mfunctionydot vdpol t y ydot zeros size y ydot 1 y 2 ydot 2 y 2 y 1 2 1 y 1 或寫為 ydot y 1 y 2 y 1 2 1 y 1 四 Matlab對應命令 ode23 ode45 調用格式 t y ode23 函數文件名 t0 tN y0 tol t y ode45 函數文件名 t0 tN y0 tol 默認精度 ode23 1e 3ode45 1e 6 說明 t0 初始時刻 tN 終點時刻y0 初值 tol 計算精度 3月15日作業(yè) 1 VanderPol方程的兩種解法 1 采用ode45命令2 Runge Kutta方法2 Duffing方程的求解 Runge Kutta方法 計算步長h 0 005 計算時間t0 0 0 tN 100 要求 寫出程序體 打印所繪圖形 圖形標題用個人的名字 Duffing方程 五 動力學系統(tǒng)的求解 1 動力學方程 2 二階方程轉成一階方程 1 令 2 其中 即 2 3 Matlab程序 主程序 ZCX t0 Y0 h N P0 w 輸入初始值 步長 迭代次數 初始激勵力 fori 1 Nt1 t0 hP 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 輸入 t0 h 時刻的外部激勵力K4 ZCX sub t0 h y0 hK3 P Y1 y0 h 6 K1 2K2 2K3 K4 t1 Y1 輸出t1 y1 nexti輸出數據或圖形 Matlab程序 子程序 ZCX sub m functionydot f t Y P M K C 輸入結構參數P1 zeros 3 1 inv M P A zeros 0 0 eye n n M 1K M 1C ydot AY P1 例題2 三自由度質量彈簧系統(tǒng) 矩陣表示 其中 動力學方程 解析解 已知參數 m1 m2 m3 1 k1 2 k2 2 K3 1 K4 2 P0 1 要求 采用四階龍格 庫塔法編程計算三個質量的響應時程 計算時間0 50 例如 4階龍格 庫塔法的結果 ode45的結果 第一個質量的位移響應時程 結果完全一致 MATLAB程序 1 4階RK方法 2 采用ode45 m chap2 ex2 1 m m chap2 ex2 1 sub m 例題3 蹦極跳系統(tǒng)的動態(tài)仿真 蹦極者系著一根彈性繩從高處的橋梁 或山崖等 向下跳 在下落的過程中 蹦極者幾乎處于失重狀態(tài) 按照牛頓運動規(guī)律 自由下落的物體由下式確定 其中 m為人體的質量 g為重力加速度 x為物體的位置 第二項和第三項表示空氣的阻力 其中位置x的基準為橋梁的基準面 即選擇橋梁作為位置的起點x 0 低于橋梁的位置為正值 高于橋梁的位置為負值 如果人體系在一個彈性常數為k的彈性繩索上 定義繩索下端的初始位置為0 則其對落體位置的影響為 空氣的阻力 整個蹦極系統(tǒng)的數學模型為 設橋梁距離地面為50m 即h2 50 蹦極者的起始位置為繩索的長度30m 即h1 30 蹦極者起始速度為0 其余的參數分別為k 20 a2 a1 1 m 70kg g 10m s2 初始條件 已知參數 初始條件變?yōu)?y0 30 0 初始位移和初始速度 t y ode45 bengji sub 0 0 01 100 y0 x1 50 y 1 x1代表蹦極者與地面之間的距離plot t x1 gridplot t y 1 grid y 1 代表位移 主程序 程序名 bengji m Matlab程序 functionydot f t y m 70 k 20 a1 1 a2 1 g 10 x y 1 x代表蹦極者的位移x dot y 2 x dot代表x的速度ifx 0ydot 0 1 k m a1 m a2 m abs x dot y 0 g elseydot 0 1 0 a1 m a2 m abs x dot y 0 g end 子程序 程序名 bengji sub m y 1 x1 結果分析 右上圖為蹦極者與地面之間的距離 從結果可看出 對于體重為70kg