Programming and Application(编程与应用)

Content(目录)

Linux

MySQL
Office

Numeric Integration Method : 4th Order Runga-Kuta

Zhanshan Dong

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²)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²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');

Post comments(留言)

由Google提供