Both interval and step oriented modes of operation are available and also modes designed to permit intermediate output within an interval oriented mode.
An outline of a typical calling program for
nag_ode_ivp_stiff_imp_sparjac (d02nj) is given below. It calls
the sparse matrix linear algebra setup function
nag_ode_ivp_stiff_sparjac_setup (d02nu),
the Backward Differentiation Formula (BDF) integrator setup
function
nag_ode_ivp_stiff_bdf (d02nv), its diagnostic counterpart
nag_ode_ivp_stiff_integ_diag (d02ny),
and
the sparse matrix linear algebra diagnostic function
nag_ode_ivp_stiff_sparjac_diag (d02nx).
.
.
.
[...] = d02nv(...);
[...] = d02nu(...);
[..., ifail] = d02nj(...);
if (ifail ~= 1 and ifail < 14)
[...] = d02nx(...);
[...] = d02ny(...);
end
.
.
.
The linear algebra setup function
nag_ode_ivp_stiff_sparjac_setup (d02nu) and
one of
the integrator setup
functions,
nag_ode_ivp_stiff_dassl (d02mv),
nag_ode_ivp_stiff_bdf (d02nv) or
nag_ode_ivp_stiff_blend (d02nw),
must be called prior to the call of
nag_ode_ivp_stiff_imp_sparjac (d02nj).
Either or both of the integrator diagnostic function
nag_ode_ivp_stiff_integ_diag (d02ny), or the
sparse matrix linear algebra diagnostic function
nag_ode_ivp_stiff_sparjac_diag (d02nx), may be called after the call to
nag_ode_ivp_stiff_imp_sparjac (d02nj).
There is also a function,
nag_ode_ivp_stiff_contin (d02nz), designed to permit you to change step size on a continuation call to
nag_ode_ivp_stiff_imp_sparjac (d02nj) without restarting the integration process.
This example solves the well-known stiff Robertson problem written as a mixed differential/algebraic system in implicit form
exploiting the fact that, from the initial conditions
$a=1.0$ and
$b=c=0.0$, we know that
$a+b+c=1$ for all time. We integrate over the range
$\left[0,10.0\right]$ with vector relative error control and scalar absolute error control (
${\mathbf{itol}}=3$) and using the BDF integrator (setup function
nag_ode_ivp_stiff_bdf (d02nv))
and a modified Newton method.
The Jacobian is evaluated, in turn, using the 'A' (Analytical) and 'F' (Full information) options.
We provide a monitor function to terminate the integration when the value of the component a falls below
$0.9$.
function d02nj_example
fprintf('d02nj example results\n\n');
global ncall tkeep ykeep jacobian;
jacobian = 'Analytic';
d02nj_subexample;
jacobian = 'Full_info';
d02nj_subexample;
fig1 = figure;
display_plot(tkeep, ykeep);
function d02nj_subexample
global ncall tkeep ykeep jacobian;
neq = int64(3);
neqmax = neq;
maxord = int64(5);
sdysav = maxord+1;
petzld = true;
const = zeros(6, 1);
tcrit = 0;
hmin = 1.0e-10;
hmax = 10;
h0 = 0;
maxstp = int64(200);
mxhnil = int64(5);
rwork = zeros(50+4*neq, 1);
[const, rwork, ifail] = d02nv(...
neqmax, sdysav, maxord, 'Newton', petzld, ...
const, tcrit, hmin, hmax, h0, maxstp, ...
mxhnil, 'Average-L2', rwork);
nwkjac = int64(100);
ia = int64([1; 3; 6; 9]);
ja = int64([1; 2;
1; 2; 3;
1; 2; 3]);
njcpvt = int64(150);
sens = 1.0e-6;
u = 0.1;
eta = 1.0e-4;
lblock = true;
[jacpvt, rwork, ifail] = d02nu(...
neq, neqmax, jacobian, nwkjac, ia, ja, ...
njcpvt, sens, u, eta, lblock, rwork);
t = 0;
tout = 10.0;
y = [1; 0; 0];
ydot = zeros(neq, 1);
itol = int64(3);
rtol = [0.0001; 0.001; 0.0001];
atol = [1e-07];
inform = zeros(23,1,'int64');
ysave = zeros(neq, sdysav);
wkjac = zeros(nwkjac, 1);
lderiv = [false; true];
itrace = int64(0);
itask = int64(1);
fprintf('Jacobian = %s\n\n',jacobian);
fprintf(' x y\n');
fprintf('%8.3f %8.5f %8.5f %8.5f\n',t,y);
if (jacobian(1)=='F')
ncall = 2;
ykeep = y;
tkeep = [t];
end
wstat = warning();
warning('OFF');
[t, tout, y, ydot, rwork, inform, ysave, wkjac, jacpvt, lderiv, ifail] ...
= d02nj(...
t, tout, y, ydot, rwork, rtol, atol, itol, inform, @resid, ...
ysave, @jac, wkjac, jacpvt, @monitr, lderiv, itask, itrace);
warning(wstat);
fprintf('%8.3f %8.5f %8.5f %8.5f\n',t,y);
[hu, h, tcur, tolsf, nst, nre, nje, nqu, nq, niter, imxer, algequ, ifail] ...
= d02ny(...
neq, neqmax, rwork, inform);
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,niter);
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);
icall = int64(0);
[liwreq, liwusd, lrwreq, lrwusd, nlu, nz, ngp, isplit, igrow, nblock]...
= d02nx(...
icall, lblock, inform);
fprintf(' sparse jacobian statistics:\n');
fprintf(' njcpvt - required = %3d, used = %3d\n',liwreq,liwusd);
fprintf(' nwkjac - required = %3d, used = %3d)\n',lrwreq,lrwusd);
fprintf(' No. of LU-decomps = %3d, No. of nonzeros = %3d\n',nlu,nz);
fprintf([' No. of function calls to form Jacobian = %3d, try ', ...
'isplit = %3d\n'], ngp, isplit);
fprintf([' Growth estimate = %d, No. of blocks on diagonal %3d\n\n'], ...
igrow, nblock);
function pdj = jac(neq, t, y, ydot, h, d, j, pdj)
pdj = zeros(neq,1);
hxd = h*d;
if j == 1
pdj(1) = -hxd*(1.0);
pdj(2) = -hxd*(0.04);
elseif j == 2
pdj(1) = -hxd*(1.0);
pdj(2) = -hxd*(-1.0e4*y(3)-6.0e7*y(2)) + 1;
pdj(3) = -hxd*(6.0e7*y(2));
elseif j == 3
pdj(1) = -hxd*(1.0);
pdj(2) = -hxd*(-1.0e4*y(2));
pdj(3) = -hxd*(0.0) + 1;
end
function [hnext, y, imon, inln, hmin, hmax] = ...
monitr(neq, neqmax, t, hlast, hnext, y, ydot, ...
ysave, r, acor, imon, hmin, hmax, nqu)
global ncall tkeep ykeep jacobian;
if (y(1) <= 0.9)
imon = int64(-2);
end
inln = int64(0);
if (jacobian(1)=='F')
tkeep(ncall,1) = t;
ykeep(:,ncall) = y;
ncall = ncall + 1;
end
function [r, ires] = resid(neq, t, y, ydot, ires)
r = zeros(neq,1);
r(1) = 0;
r(2) = -ydot(2);
r(3) = -ydot(3);
if ires == 1
r(1) = y(1) + y(2) + y(3) - 1 + r(1);
r(2) = 0.04*y(1) - 1.0e4*y(2)*y(3) - 3.0e7*y(2)*y(2) + r(2);
r(3) = 3.0e7*y(2)*y(2) + r(3);
end
function display_plot(x, y)
hline1 = semilogx(x, y(1,:));
hold on
[haxes, hline2, hline3] = plotyy(x, y(3,:),x, y(2,:), ...
'semilogx', 'semilogx');
set(haxes(1), 'YLim', [-0.05 1.2]);
set(haxes(1), 'XMinorTick', 'on', 'YMinorTick', 'on');
set(haxes(1), 'YTick', [0.0:0.2:1.2]);
set(haxes(2), 'YLim', [0.0 4e-5]);
set(haxes(2), 'YMinorTick', 'on');
set(haxes(2), 'YTick', [5e-6:5e-6:3.5e-5]);
for iaxis = 1:2
set(haxes(iaxis), 'XLim', [0 12]);
set(haxes(iaxis), 'XTick', [0.0001 0.001 0.01 0.1 1 10]);
end
set(gca, 'box', 'off');
ht = title({'Stiff Robertson Problem (Implicit DAE):',...
'BDF with Newton Iterations'});
set(ht,'Position',[0.02,1.15,0.0]);
xlabel('x');
ylabel(haxes(1),'Solution (a,c)');
ylabel(haxes(2),'Solution (b)');
legend('a','c','b','Location','Best');
set(hline1, 'Marker', '+','Linestyle','-','Color','red');
set(hline2, 'Marker', '*','Linestyle',':','Color','green');
set(hline3, 'Marker', 'x','Linestyle','--','Color','blue');