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_exp_bandjac (d02nc) is given below. It calls the banded matrix linear algebra setup function
nag_ode_ivp_stiff_bandjac_setup (d02nt), the Backward Differentiation Formula (BDF) integrator setup function
nag_ode_ivp_stiff_bdf (d02nv), and its diagnostic counterpart
nag_ode_ivp_stiff_integ_diag (d02ny).
.
.
.
[...] = d02nv(...);
[...] = d02nt(...);
[..., ifail] = d02nc(...);
if (ifail ~= 1 and ifail < 14)
[...] = d02ny(...);
end
.
.
.
The linear algebra setup function
nag_ode_ivp_stiff_bandjac_setup (d02nt) and one of the integrator setup functions,
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_exp_bandjac (d02nc). The integrator diagnostic function
nag_ode_ivp_stiff_integ_diag (d02ny) may be called after the call to
nag_ode_ivp_stiff_exp_bandjac (d02nc). 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_exp_bandjac (d02nc) without restarting the integration process.
The cost of computing a solution depends critically on the size of the differential system and to a lesser extent on the degree of stiffness of the problem. For nag_ode_ivp_stiff_exp_bandjac (d02nc) the cost is proportional to ${\mathbf{neq}}\times {\left({\mathbf{ml}}+{\mathbf{mu}}+1\right)}^{2}$ though for problems which are only mildly nonlinear the cost may be dominated by factors proportional to ${\mathbf{neq}}\times \left({\mathbf{ml}}+{\mathbf{mu}}+1\right)$ except for very large problems.
function d02nc_example
fprintf('d02nc example results\n\n');
global ncall ykeep tkeep
neq = int64(3);
neqmax = int64(neq);
ml = int64(1);
mu = int64(2);
nwkjac = int64(neqmax*(2*ml+mu+1));
njcpvt = int64(neqmax);
maxord = int64(11);
sdysav = int64(maxord+3);
maxstp = int64(200);
mxhnil = int64(5);
h0 = 0;
hmax = 10;
hmin = 1.0e-10;
tcrit = 0;
petzld = false;
const = zeros(6, 1);
rwork = zeros(50+4*neqmax, 1);
[const, rwork, ifail] = d02nw(...
neqmax, sdysav, maxord, const, tcrit, hmin, ...
hmax, h0, maxstp, mxhnil, 'Average-L2', rwork);
[rwork, ifail] = d02nt(...
neq, neqmax, 'Analytic', ml, mu, nwkjac, njcpvt, rwork);
inform(1:32) = int64(0);
ysave = zeros(neq, sdysav);
wkjac = zeros(nwkjac, 1);
jacpvt(1:njcpvt) = int64(0);
t = 0.0;
tout = 5.0;
itask = int64(1);
itrace = int64(0);
y = [1; 0; 0];
itol = int64(2);
rtol = [0.0001];
atol = [1e-07; 1e-08; 1e-07];
ncall = 1;
tkeep(1) = t;
ykeep(1:3,1) = y(1:3);
fprintf(' x y_1 y_2 y_3\n');
fprintf('%8.3f %8.5f %8.5f %8.5f\n', t, y);
[t, y, ydot, rwork, inform, ysave, wkjac, jacpvt, ifail] = ...
d02nc(t, tout, y, rwork, rtol, atol, itol, inform, @fcn, ysave, @jac, ...
wkjac, jacpvt, @monitr, itask, itrace);
fprintf('%8.3f %8.5f %8.5f %8.5f\n', t, y);
h = 0.7;
[rwork, ifail] = d02nz(neqmax, tcrit, h, hmin, hmax, maxstp, mxhnil, rwork);
tout = 10;
[t, y, ydot, rwork, inform, ysave, wkjac, jacpvt, ifail] = ...
d02nc(t, tout, y, rwork, rtol, atol, itol, inform, @fcn, ysave, @jac, ...
wkjac, jacpvt, @monitr, itask, itrace);
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\n',imxer);
fig1 = figure;
display_plot(tkeep(1:ncall),ykeep(1:3,1:ncall));
function [f, ires] = fcn(neq, t, y, ires)
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, ml, mu, pIn)
p = zeros(ml+mu+1, neq);
hxd = h*d;
p(1,1) = -hxd*(-0.04d0) + 1;
p(2,1) = -hxd*(1.0d4*y(3));
p(3,1) = -hxd*(1.0d4*y(2));
p(1,2) = -hxd*(0.04d0);
p(2,2) = -hxd*(-1.0d4*y(3)-6.0d7*y(2)) + 1;
p(3,2) = -hxd*(-1.0d4*y(2));
p(1,3) = -hxd*(6.0d7*y(2));
p(2,3) = -hxd*(0.0d0) + 1;
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
tend = 10.0;
if (imon == 1 && t>0.001 && t<tend)
ncall = ncall + 1;
tkeep(ncall) = t;
ykeep(1:3,ncall) = y(1:3);
end
inln = int64(0);
function display_plot(tkeep,ykeep)
hline3 = plot(tkeep,ykeep(3,:));
hold on;
[haxes, hline1, hline2] = plotyy(tkeep,ykeep(1,:), tkeep,ykeep(2,:));
set(haxes(1), 'YLim', [0.0 1.1]);
set(haxes(1), 'XMinorTick', 'on', 'YMinorTick', 'on');
set(haxes(1), 'YTick', [0.0 0.2 0.4 0.6 0.8 1.0]);
set(haxes(2), 'YLim', [0.0 4e-5]);
set(haxes(2), 'YMinorTick', 'on');
set(haxes(2), 'YTick', [0.5e-5:0.5e-5:3.5e-5]);
set(haxes(1), 'XLim', [0 9]);
set(haxes(2), 'XLim', [0 9]);
set(gca, 'box', 'off');
ht = title({'Stiff Robertson Problem:','BLEND integrator, Full Jacobian'});
set(ht,'Position',[4.5,0.99,0.0]);
xlabel('x');
ylabel(haxes(1),'Solution (a,c)');
ylabel(haxes(2),'Solution (b)');
legend('c','a','b','Location','West');
set(hline1, 'Marker', '+','Linestyle','-','Color','red');
set(hline2, 'Marker', '*','Linestyle',':','Color','blue');
set(hline3, 'Marker', 'x','Linestyle','--','Color','green');
hold off;