微分方程數(shù)值解法matlab(四階龍格-庫塔法)

微分方程數(shù)值解法matlab(四階龍格—庫塔法)第一頁，共36頁。

微分方程的數(shù)值解法四階龍格—庫塔法（TheFourth-OrderRunge－KuttaMethod）第二頁，共36頁。常微分方程（Ordinarydifferentialequations,ODE)初值問題---給出初始值邊值問題---給出邊界條件與初值常微分方程解算有關(guān)的指令ode23ode45

ode113ode23tode15sode23sode23tb第三頁，共36頁。一.解ODE的基本機理:2.把高階方程轉(zhuǎn)換成一階微分方程組1.列出微分方程初始條件令(2.1)(2.2)(2.3)第四頁，共36頁。例：著名的VanderPol方程

令降為一階初始條件第五頁，共36頁。3.根據(jù)式(2.2)編寫計算導(dǎo)數(shù)的M函數(shù)文件-ODE文件把t,Y作為輸入宗量，把作為輸出宗量%Mfunctionfilename:dYdt.m

function

Yd=f(t,Y)

Yd=f(t,Y)

的展開式例VanderPol方程%Mfunctionfilename:dYdt.m

function

Yd=f(t,Y)

Yd=zeros(size(Y));第六頁，共36頁。4.使編寫好的ODE函數(shù)文件和初值供微分方程解算指令（solver）調(diào)用Solver解算指令的使用格式[t,Y]=solver(‘ODE函數(shù)文件名’,t0,tN,Y0,tol);ode45輸出宗量形式說明：t0：初始時刻；tN：終點時刻Y0：初值；tol：計算精度第七頁，共36頁。例題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命令第八頁，共36頁。VanderPol方程

%子程序

（程序名：dYdt.m）

functionYdot=dYdt(t,Y)Ydot=[Y(2);-Y(2)*(Y(1)^2-1)-Y(1)];或?qū)憺閒unctionYdot=dYdt(t,Y)Ydot=zeros(size(Y));Ydot(1)=Y(2);Ydot(2)=-Y(2)*(Y(1).^2-1)-Y(1)];第九頁，共36頁。第十頁，共36頁。解法指令解題類型特點適合場合ode45非剛性采用4、5階Runge－Kutta法大多數(shù)場合的首選算法ode23非剛性采用Adams算法較低精度（10－3）場合ode113非剛性多步法；采用Adams算法；高低精度均可（10－3～10－6）ode45計算時間太長時取代ode45ode23t適度剛性采用梯形法則算法適度剛性ode15s剛性多步法；采用2階Rosenbrock算式，精度中等當(dāng)ode45失敗時使用；或存在質(zhì)量矩陣時ode23s剛性一步法；采用2階Rosenbrock算式，低精度低精度時，比ode15s有效；或存在質(zhì)量矩陣時ode23tb剛性采用梯形法則－反向數(shù)值微分兩階段算法，低精度低精度時，比ode15s有效；或存在質(zhì)量矩陣時各種solver解算指令的特點第十一頁，共36頁。二.四階Runge-Kutta法對I=[a,b]作分割步長第十二頁，共36頁。初值問題的數(shù)值解法分為兩大類單步法-Runge-Kutta方法多步法-Admas方法計算的近似值時只用到，是自開始方法Runge-Kutta法是常微分方程的一種經(jīng)典解法MATLAB對應(yīng)命令：ode45第十三頁，共36頁。四階Runge-Kutta公式第十四頁，共36頁。四階Runge-Kutta法計算流程圖開始Nextifori=1:N

Plot初始條件：；

積分步長：迭代次數(shù)：輸出結(jié)果子程序計算End第十五頁，共36頁。三.Runge-Kutta法解VanderPol方程的Matlab程序結(jié)構(gòu)

主程序：RK_vanderpol.m

子程序：RK_sub.m（函數(shù)文件）第十六頁，共36頁。解法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第十七頁，共36頁。第十八頁，共36頁。子程序：RK_sub.m

functionydot=vdpol(t,y)

ydot=zeros(size(y));ydot(1)=y(2);

ydot(2)=-y(2)*(y(1)^2-1)-y(1);

或?qū)憺椋?/p>

ydot=[y(1);-y(2)*(y(1)^2-1)-y(1)];第十九頁，共36頁。四.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)精度：ode23——1e-3ode45——1e-6說明：t0：初始時刻；tN：終點時刻y0：初值；tol：計算精度第二十頁，共36頁。3月15日作業(yè):1.VanderPol方程的兩種解法:1)采用ode45命令2)Runge-Kutta方法2.Duffing方程的求解（Runge-Kutta方法，計算步長h=0.005，計算時間t0=0.0,tN=100）要求：寫出程序體，打印所繪圖形，圖形標(biāo)題用個人的名字。Duffing方程第二十一頁，共36頁。五.動力學(xué)系統(tǒng)的求解1.動力學(xué)方程2.二階方程轉(zhuǎn)成一階方程(1)令：(2)第二十二頁，共36頁。其中：即：(2)第二十三頁，共36頁。3.Matlab程序(主程序：ZCX)

t0;Y0;h;N;P0,w;%輸入初始值、步長、迭代次數(shù)、初始激勵力；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輸出數(shù)據(jù)或圖形第二十四頁，共36頁。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第二十五頁，共36頁。例題2：三自由度質(zhì)量彈簧系統(tǒng)m1m2m3k1k2k3x1x2x3k4P0sin(wt)第二十六頁，共36頁。矩陣表示其中：第二十七頁，共36頁。動力學(xué)方程：解析解：已知參數(shù)：m1=m2=m3=1,k1=2,k2=2,K3=1,K4=2,P0=1,要求：采用四階龍格－庫塔法編程計算三個質(zhì)量的響應(yīng)時程.計算時間0~50例如：第二十八頁，共36頁。4階龍格－庫塔法的結(jié)果ode45的結(jié)果第一個質(zhì)量的位移響應(yīng)時程結(jié)果完全一致MATLAB程序（1）4階RK方法：

（2）采用ode45：m_chap2_ex2_1.m，m_chap2_ex2_1_sub.m第二十九頁，共36頁。例題3：蹦極跳系統(tǒng)的動態(tài)仿真蹦極者系著一根彈性繩從高處的橋梁（或山崖等）向下跳。在下落的過程中，蹦極者幾乎處于失重狀態(tài)。按照牛頓運動規(guī)律，自由下落的物體由下式確定：其中，m為人體的質(zhì)量，g為重力加速度，x為物體的位置，第二項和第三項表示空氣的阻力。其中位置x的基準(zhǔn)為橋梁的基準(zhǔn)面（即選擇橋梁作為位置的起點x＝0），低于橋梁的位置為正值，高于橋梁的位置為負值。如果人體系在一個彈性常數(shù)為k的彈性繩索上，定義繩索下端的初始位置為0，則其對落體位置的影響為：地面x橋梁基準(zhǔn)面0梯子h2h1空氣的阻力第三十頁，共36頁。整個蹦極系統(tǒng)的數(shù)學(xué)模型為：設(shè)橋梁距離地面為50m，即h2=50，蹦極者的起始位置為繩索的長度30m，即h1=30，蹦極者起始速度為0，其余的參數(shù)分別為k＝20，a2＝a1＝1；m＝70kg，g＝10m/s2。地面x橋梁基準(zhǔn)面0梯子h2h1初始條件：已知參數(shù)：第三十一頁，共36頁。令：初始條件變?yōu)椋旱谌?，?6頁。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程序第三十三頁，共36頁。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;-