太阳帆数据 Sunfine Data

Numeric Integration Method : 4th Order Runga-Kuta

By Zhanshan Dong on August 26th, 2009

ODE methods in MATLAB is so slow. I could not stant them when I solved a big problem. First I used Euler method to do numeric intergration. My advisor just laughed at me about it. He wanted me use Runga-Kuta method. I know ODE solvers in MATLAB already use Runga-Kuta methods. I just think it is slow. I could not use it. But I tried to write a simple one do not include the gabage code. I rewrote the Runga-Kuta integration method in Matlab. It works just fine. If the delta t that use coose in oderk4 is proper, the result just like the output of ODE45 which use the fourth oder Runga-Kuta method and with a procedure to find the right step size.

function [T,Y] = oderk4(funname, tspan, y)
dt = 0.1;
T  = [];
Y  = [];
T  = [T tspan(1)];
Y  = [Y y];

for t = tspan(1):dt:tspan(2)-dt
    dt2 = dt/2.0;
    dt6 = dt/6.0;
    xt  = t + dt2;

    dydt = feval(funname, t, y);
    tmpy = y + dt2 * dydt;

    dy4  = feval(funname, xt, tmpy);
    tmpy = y + dt2 * dy4;

    dy3  = feval(funname, xt, tmpy);
    tmpy = y + dt * dy3;
    dy3  = dy3 + dy4;

    dy4  = feval(funname, t + dt, tmpy);
    y    = y + dt6 * (dydt + dy4 + 2.0*dy3);

    T = [T t+dt];
    Y = [Y y];
end

van der Pol equation:

 y1"-mu(1-y1^2)y1'+y1=0 

where mu >0 is a scalar parameter.
The equation can be reeriten as a first-order differential equations:

y1'=y2
y2'=mu(1-y1^2)y2 - y1

The ODE file is as follows.

function dy=vdp1(t,y)
dy=[y(2);(1-y(1)^2)*y(2)-y(1)];

To solve this equation by Runga-Kuta procedure by using the following MATLAB command:

[T, Y]=oderk4(@vdp1,[0,30],[2;0]);

Using the following command to plot the result.

plot(T,Y(1,:),'-',T,Y(2,:),'--');
title('Solution of van der Pol Equation, mu=1');
xlabel('time t');
ylabel('solution y');
legend('y1','y2');

The result is as same as the ODE45 in MATLAB. Click the following link to show the figure.