oderk6.m

FEM tools for caculation of nonlinear problems

function [time, y_sol] = odeRK6( dydt, time, y0, u, epps )

global LK

if nargin < 4
    error(' Not enough input arguments.  See odeRK6 ');
end

 warning off;
 pts    = length(time);
 n      = length(y0);
 y_sol  = zeros(pts,3*LK); 

 y0 = y0(:);  
 y_sol(1,1:n) = y0';
 y_old = y0;


 a1=0;  a2 = 1/5;    a3 = 3/10;  a4 = 3/5;   a5 = 1;     a6 = 7/8; 
 b21= 1/5;
 b31= 3/40;   b32= 9/40;
 b41= 3/10;   b42=-9/10;  b43= 6/5;
 b51= -11/54; b52=5/2;    b53=-70/27; b54= 35/27;
 b61= 1631/55296; b62= 175/512; b63=575/13824; b64=44275/110592; b65= 253/4096;
 
 for p = 2:pts
     passit = 0;
     t      = time(p);
     dt     = t - time(p-1);
     t0     = t - dt; TT0    = t0;
     YOLD   = y_old;
     hmin = dt/100000;
     
     DTT = dt;
     stN = 1;                       
     goOUT = 0;

     u0   = u(:,p-1);              
     dudt = (u(:,p)-u0)/dt;        
     
     while passit==0   

           
         k1 = feval( dydt, t0 + DTT*a1, y_old, u0+dudt*(t0+DTT*a1-t0) );
         k2 = feval( dydt, t0 + DTT*a2, y_old + k1*DTT*b21, u0+dudt*(t0+DTT*a2-t0) );
         k3 = feval( dydt, t0 + DTT*a3, y_old + k1*DTT*b31 + k2*DTT*b32, u0+dudt*(t0+DTT*a3-t0) );
         k4 = feval( dydt, t0 + DTT*a4, y_old + k1*DTT*b41 + k2*DTT*b42 + k3*DTT*b43, u0+dudt*(t0+DTT*a4-t0) );
         k5 = feval( dydt, t0 + DTT*a5, y_old + k1*DTT*b51 + k2*DTT*b52 + k3*DTT*b53 + k4*DTT*b54, u0+dudt*(t0+DTT*a5-t0)  );
         k6 = feval( dydt, t0 + DTT*a6, y_old + k1*DTT*b61 + k2*DTT*b62 + k3*DTT*b63 + k4*DTT*b64 + k5*DTT*b65, u0+dudt*(t0+DTT*a6-t0) );


           y_n5   = y_old + ( k1*37/378 + k3*250/621 + k4*125/594 + k6*512/1771 )*DTT;
           y_n4   = y_old + ( k1*2825/27648 + k3*18575/48384 + k4*13525/55296 + k5*277/14336 + k6*1/4 )*DTT;
           EE = y_n5 - y_n4;

           wish_eps = y_n5.*epps;  
           if EE == 0 EE=wish_eps; end; 

           stN = stN - 1;          

           y_old = y_n5;
           t0 = t0 + DTT;

           if goOUT==1
               if stN==0
                   disp(['   ---  step OK  ---             ' num2str(p) ]); 
                   passit = 1; 
               end;
           else
               if stN==0
                   step_INC = (max(abs(wish_eps./EE) ))^0.2 ; 
                   if step_INC == 0   step_INC=1;   end;
                   if isnan(step_INC) step_INC=1; end
                   step_INC = 1*step_INC; 
                   DTT=1.00*DTT*step_INC; t0 = TT0; y_old = YOLD; stN=ceil(dt/DTT); DTT=dt/stN;
                   if DTT<hmin 
                       error(['integration step is smaller than min -- ' num2str(hmin) ]);
                   end;
                   if isnan(stN)>=1 
                       error('  integration is not possible -- NaN --  may be too stiff');
                   end
                   if stN==0 passit = 1; end
                   disp(['   steps: ' num2str(stN)]); 
                   goOUT = 1;   
               end;
           end;


       end;

        y_sol(p,1:n) = y_old';
   end
warning on;