MATLAB龍格庫塔方法解微分方程

1、龍格-庫塔方法是一種經(jīng)典方法，具有很高的精度，它間接的利用了泰勒級數(shù)展開，避免了高階偏導(dǎo)數(shù)的計(jì)算。此處以最為經(jīng)典的四級四階龍格-庫塔方法為例，計(jì)算格式如下 1龍格-庫塔法解一階ODE對于形如 的一階ODE初值問題，可以直接套用公式，如今可以借助計(jì)算機(jī)方便的進(jìn)行計(jì)算，下面給出一個(gè)實(shí)例 取步長h=0.1，此處由數(shù)學(xué)知識可得該方程的精確解為 。在這里利用MATLAB編程，計(jì)算數(shù)值解并與精確解相比，代碼如下：(1)寫出微分方程，便于調(diào)用和修改function val = odefun( x,y ) val = y-2*x/y; end(2)編寫runge-kutta方法的函數(shù)代碼function y

2、= runge_kutta( h,x0,y0 )k1 = odefun(x0,y0);k2 = odefun(x0+h/2,y0+h/2*k1);k3 = odefun(x0+h/2,y0+h/2*k2);k4 = odefun(x0+h,y0+h*k3); y = y0+h*(k1+2*k2+2*k3+k4)/6; end(3)編寫主函數(shù)解微分方程,并觀察數(shù)值解與精確解的差異clear allh = 0.1;x0 = 0;y0 = 1;x = 0.1:h:1;y(1) = runge_kutta(h,x0,y0);for k=1:length(x) x(k) = x0+k*h; y(k+1)

3、 = runge_kutta(h,x(k),y(k);endz = sqrt(1+2*x);plot(x,y,*);hold onplot(x,z,'r');結(jié)果如下圖，數(shù)值解與解析解高度一致2龍格-庫塔法解高階ODE對于高階ODE來說，通用的方法是將高階方程通過引入新的變量降階為一階方程組，此處仍以一個(gè)實(shí)例進(jìn)行說明。 這是一個(gè)二階ODE,描述的是一個(gè)物體的有阻尼振動(dòng)情況。初始條件為 ，將方程降階，引入一個(gè)向量型變量Y故有 記則至此，二階方程降階為一階方程組。值得注意的是此時(shí)再用龍格-庫塔法進(jìn)行求解時(shí)，代入的將是一個(gè)Y向量。同樣利用MATLAB進(jìn)行計(jì)算，步長h=0.05，時(shí)間周

4、期為0,20.(1) 編寫ODE函數(shù)function Y = odefun1( ,Y0 )% 此處Y0為一個(gè)列向量，因?yàn)闀r(shí)間t未顯含在一階方程組中% 所以ode函數(shù)的第一個(gè)參數(shù)為空，要根據(jù)具體情況而定。Y = Y0(2); (2000-200*Y0(2)-750*Y0(1)/500;end(2) 編寫runge-kutta函數(shù)function Y = rkfa( h,t0,Y0 )k1 = odefun1(t0,Y0);k2 = odefun1(t0+h/2,Y0+h/2*k1);k3 = odefun1(t0+h/2,Y0+h/2*k2);k4 = odefun1(t0+h,Y0+h*k3);Y = Y0+h*(k1+2*k2+2*k3+k4)/6;end(3) 編寫主函數(shù)clear allh = 0.05;t = 0.05:h:20;t0 = 0;Y0 = 0; 0;%初值Y = cell(1,length(t);Y1 = rkfa( h,t0,Y0 );z = zeros(2,length(t);for k=1:length(t) Yk+1 = rkfa( h,t0,Yk); z(1,k) = Yk(1);z(2,k) = Yk(2