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_fulljac (d02ng) is given below. It calls the full matrix linear algebra setup function
nag_ode_ivp_stiff_fulljac_setup (d02ns), 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(...);
[...] = d02ns(...);
[..., ifail] = d02ng(...);
if (ifail ~= 1 and ifail < 14)
[...] = d02ny(...);
end
.
.
.
The linear algebra setup function
nag_ode_ivp_stiff_fulljac_setup (d02ns) 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_fulljac (d02ng). The integrator diagnostic function
nag_ode_ivp_stiff_integ_diag (d02ny) may be called after the call to
nag_ode_ivp_stiff_imp_fulljac (d02ng). 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_fulljac (d02ng) without restarting the integration process.
The accuracy of the numerical solution may be controlled by a careful choice of the arguments
rtol and
atol, and to a much lesser extent by the choice of norm. You are advised to use scalar error control unless the components of the solution are expected to be poorly scaled. For the type of decaying solution typical of many stiff problems, relative error control with a small absolute error threshold will be most appropriate (that is, you are advised to choose
${\mathbf{itol}}=1$ with
${\mathbf{atol}}\left(1\right)$ small but positive).
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_imp_fulljac (d02ng) the cost is proportional to ${{\mathbf{neq}}}^{3}$, though for problems which are only mildly nonlinear the cost may be dominated by factors proportional to ${{\mathbf{neq}}}^{2}$ except for very large problems.
function d02ng_example
fprintf('d02ng example results\n\n');
istat = 1;
neq = int64(3);
neqmax = neq;
maxord = int64(5);
sdysav = maxord+1;
petzld = false;
tcrit = 0;
hmin = 1.0e-10;
hmax = 10;
h0 = 0;
maxstp = int64(200);
mxhnil = int64(5);
const = zeros(6, 1);
rwork = zeros(50+4*neq, 1);
[const, rwork, ifail] = ...
d02nv( ...
neqmax, sdysav, maxord, 'Functional', petzld, ...
const, tcrit, hmin, hmax, h0, maxstp, mxhnil, 'Average-L2', rwork);
nwkjac = int64(neq*(neq+1));
[rwork, ifail] = d02ns( ...
neq, neqmax, 'Numerical', nwkjac, rwork);
sol = zeros(neq, 1);
wkjac = zeros(nwkjac, 1);
ysave = zeros(neq, sdysav);
ydot = zeros(neq, 1);
inform(1:23) = int64(0);
algequ = false(neq, 1);
t = 0;
tout = 0.1;
itask = int64(2);
itrace = int64(0);
y = [1; 0; 0];
lderiv = [false; false];
itol = int64(4);
rtol = [0.0001; 0.001; 0.0001];
atol = [1e-07; 1e-08; 1e-07];
fprintf(' x y_1 y_2 y_3\n');
fprintf('%8.3f %8.4f %8.6f %8.4f\n',t,y);
ncall = 1;
tkeep = t;
ykeep = y;
tstep = 0.02;
iout = 1;
xout = iout*tstep;
while t < tout
[t, tout, y, ydot, rwork, inform, ysave, wkjac, lderiv, ifail] = ...
d02ng( ...
t, tout, y, ydot, rwork, rtol, atol, itol, inform, @resid, ...
ysave, 'd02ngz', wkjac, 'd02nby', lderiv, itask, itrace);
if (t<tstep)
ncall = ncall+1;
ykeep(:,ncall) = y;
tkeep(ncall) = t;
elseif (xout <= t)
[hu, h, tcur, tolsf, nst, nre, nje, nqu, nq, niter, imxer, algequ, ...
ifail] = d02ny( ...
neq, neqmax, rwork, inform);
[sol, ifail] = d02xj( ...
xout, neq, ysave, neq, tcur, ...
nqu, hu, h, 'sdysav', sdysav);
ncall = ncall + 1;
ykeep(:,ncall) = sol;
tkeep(ncall) = xout;
fprintf('%8.3f %8.4f %8.6f %8.4f\n',xout,sol);
iout = iout + 1;
xout = iout*tstep;
end
end
if istat == 1
[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);
end
fig1 = figure;
display_plot(tkeep, ykeep)
function [r, ires] = resid(neq, t, y, ydot, ires)
r(1:3) = -ydot(1:3);
if ires == 1
r(1) = -0.04*y(1) + 1.0e4*y(2)*y(3) + 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)
hline3 = plot(x, y(3,:));
hold on;
[haxes, hline1, hline2] = plotyy(x, 100*y(2,:), x, y(1,:));
set(haxes(1), 'YLim', [0.0 0.0045]);
set(haxes(1), 'XMinorTick', 'on', 'YMinorTick', 'on');
set(haxes(1), 'YTick', [0:0.001:0.004]);
set(haxes(2), 'YLim', [0.995 1.006]);
set(haxes(2), 'YMinorTick', 'on');
set(haxes(2), 'YTick', [0.995:0.002:1.005]);
set(haxes(1), 'XLim', [-0.005 0.1]);
set(haxes(2), 'XLim', [-0.005 0.1]);
set(gca, 'box', 'off');
ht = title({'Stiff Robertson Problem:','BDF and Functional Iteration'});
set(ht,'Position',[0.05,0.004,0.0]);
xlabel('x');
ylabel(haxes(1), 'Solution (100*b,c)');
ylabel(haxes(2), 'Solution (a)');
legend('c', '100*b', 'a', 'Location', 'Best');
set(hline1,'Linewidth',0.5,'Marker','+','LineStyle','-','Color','red');
set(hline2,'Linewidth',0.5,'Marker','*','LineStyle',':','Color','blue');
set(hline3,'Linewidth',0.5,'Marker','x','LineStyle','--','Color','green');
hold off;