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

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