微分方程的數值解法 四階龍格—庫塔法 (The Fourth-Order Runge－Kutta Method) 常微分方程(Ordinary differential equations, ODE) 初值問題---給出初始值 邊值問題---給出邊界條件 一.解ODE的基本機理: 3. 根據式(2.2)編寫計算導數的M函數文件-ODE文件 例題1：著名的Van der Pol方程 四階Runge-Kutta公式 三. Runge-Kutta 法解Van der Pol 方程的Matlab 程序結構主程序：RK_vanderpol.m 子程序：RK_sub.m(函數文件) 解法2：采用Runge_Kutta法編程計算 主程序：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); 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 四. Matlab對應命令：ode23,ode45 五. 動力學系統的求解 1. 動力學方程 3. Matlab 程序(主程序：ZCX) t0;Y0;h;N;P0,w; %輸入初始值、步長、迭代次數、初始激勵力； 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= %輸入 (t0+h) 時刻的外部激勵力 K4 = ZCX_sub (t0 + h, y0 + hK3, P) Y1 = y0 + (h/6) (K1 + 2K2 + 2K3 + K4) t1, Y1 (輸出 t1, y1) next i 輸出數據或圖形 Matlab 程序(子程序：ZCX_sub.m) 例題2：三自由度質量彈簧系統 * * 與初值常微分方程解算有關的指令 ode23 ode45 ode113 ode23t ode15s ode23s ode23tb 2. 把高階方程轉換成一階微分方程組 1. 列出微分方程 初始條件 令 (2.1) (2.2) (2.3) 例：著名的Van der Pol方程 令 降為一階 初始條件 把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函數文件和初值 供微分方程解算指令(solver)調用 Solver解算指令的使用格式 [t, Y]=solver (‘ODE函數文件名’, t0, tN, Y0, tol); ode45 輸出宗量形式 說明： t0：初始時刻；tN：終點時刻Y0：初值； tol：計算精度 % 主程序

總結

以上是生活随笔為你收集整理的四阶五级matlab,微分方程数值解法matlab(四阶龙格—库塔法).ppt的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。