ode4m.m

来自「结构力学中的有限元例子,包含了7个分类文件夹」· M 代码 · 共 29 行

function [time, y_sol] = ode4m( dydt, time, y0 )
% Henri Gavin, Civil and Environmental Engineering, Duke Univsersity, Jan. 2001
% WH Press, et al, Numerical Recipes in C, 1992, Section 16.1
 pts = length(time);
 n      = length(y0);
 y_sol  = zeros(pts,n);
 
 y0 = y0(:);  
 y_sol(1,:) = y0';
 y_old = y0;
 
 for p = 2:pts

     t      = time(p); 
     dt     = t - time(p-1);
     dt2    = dt/2;

     k1  = feval( dydt, t - dt,  y_old );
     k2  = feval( dydt, t - dt2, y_old + k1*dt2 );
     k3  = feval( dydt, t - dt2, y_old + k2*dt2 );
     k4  = feval( dydt, t,       y_old + k3*dt );

     y_old = y_old + ( k1 + 2*(k2 + k3) + k4 )*dt/6; 

     y_sol(p,:) = y_old';

 end