编程技术、软件应用与系统模拟

(Programming, Applicaiton and Simulation)


		本站目录
		首页 ASP/Access/IIS DELPHI/PASCAL PASCAL高级编程 C语言编程实例 WORD Excel MATLAB MINITAB讲座 Windows DOS SAS 生物系统模拟土壤水分剖析器其他


		镜像站点
		主站北美镜象站欧洲镜象站(1) 欧洲镜象站(2)

座右铭:
只做有益人类的事
不做有害人类的事

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

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