NAG Toolbox: nag_ode_ivp_stiff_integ_diag (d02ny)

function d02ny_example


fprintf('d02ny example results\n\n');

% Initialize integration method setup variables and arrays.
neq    = int64(3);
neqmax = int64(neq);
nwkjac = int64(neqmax*(neqmax + 1));
maxord = int64(5);
sdysav = int64(maxord+1);
maxstp = int64(200);
mxhnil = int64(5);

h0    = 0;
hmax  = 10;
hmin  = 1.0e-10;
tcrit = 20;
petzld = false;

const  = zeros(6, 1);
rwork  = zeros(50+4*neqmax, 1);

[const, rwork, ifail] = d02nv(neqmax, sdysav, maxord, 'Newton', petzld, ...
                              const, tcrit, hmin, hmax, h0, maxstp, ...
                              mxhnil, 'Average-L2', rwork);

% Setup for numerical Jacbian using d02ns
nwkjac = int64(neq*(neq+1));
[rwork, ifail] = d02ns( ...
			neq, neq, 'Numerical', nwkjac, rwork);

% Initialize variables and arrays for integration
t      = 0;
tout   = 10;
y      = [1; 0; 0];
rtol   = [0.0001];
atol   = [1e-07];
itol   = int64(1);
inform = zeros(23, 1, 'int64');
ysave  = zeros(neq, sdysav);
wkjac  = zeros(nwkjac, 1);
itask  = int64(4);
itrace = int64(0);

% Integrate ODE from t=0 to t=tout, no monitoring, using d02nb
[t, y, ydot, rwork, inform, ysave, wkjac, ifail] = ...
  d02nb( ...
	 t, tout, y, rwork, rtol, atol, itol, inform, @fcn, ysave, ...
	 @jac, wkjac, 'd02nby', itask, itrace);

% Get diagnostics
[hu, h, tcur, tolsf, nst, nre, nje, nqu, nq, nit, imxer, algequ, ifail] = ...
  d02ny( ...
         neq, neqmax, rwork, inform);

% Get solution at tout
[y, ifail] = d02mz( ...
                    tout, neq, neq, ysave, rwork);

% Print solution and diagnostics
fprintf(' At t = %8.3f, y(1:3) = %5.1f, %5.1f, %5.1f\n', t, y);
fprintf('\nDiagnostic information\n integration status:\n');
fprintf('   last and next step sizes = %8.5f, %8.5f\n',hu, h);
fprintf('   integration stopped at x = %8.5f\n',tcur);
fprintf(' algorithm statistics:\n');
fprintf('   number of time-steps and Newton iterations  = %5d %5d\n',nst,nit);
fprintf('   number of residual and jacobian evaluations = %5d %5d\n',nre,nje);
fprintf('   order of method last used and next to use   = %5d %5d\n',nqu,nq);
fprintf('   component with largest error                = %5d\n',imxer);



function [f, ires] = fcn(neq, t, y, ires)
  % Evaluate derivative vector.
  f = zeros(3,1);
  f(1) = -0.04d0*y(1) + 1.0d4*y(2)*y(3);
  f(2) = 0.04d0*y(1) - 1.0d4*y(2)*y(3) - 3.0d7*y(2)*y(2);
  f(3) = 3.0d7*y(2)*y(2);

function p = jac(neq, t, y, h, d, p)
  % Evaluate the Jacobian.
  p = zeros(neq, neq);
  hxd = h*d;
  p(1,1) = 1.0d0 - hxd*(-0.04d0);
  p(1,2) = -hxd*(1.0d4*y(3));
  p(1,3) = -hxd*(1.0d4*y(2));
  p(2,1) = -hxd*(0.04d0);
  p(2,2) = 1.0d0 - hxd*(-1.0d4*y(3)-6.0d7*y(2));
  p(2,3) = -hxd*(-1.0d4*y(2));
  p(3,2) = -hxd*(6.0d7*y(2));
  p(3,3) = 1.0d0 - hxd*(0.0d0);

d02ny example results

 At t =   10.000, y(1:3) =   0.8,   0.0,   0.2

Diagnostic information
 integration status:
   last and next step sizes =  0.90178,  0.90178
   integration stopped at x = 10.76621
 algorithm statistics:
   number of time-steps and Newton iterations  =    55    78
   number of residual and jacobian evaluations =   128    16
   order of method last used and next to use   =     4     4
   component with largest error                =     3