Note: the interface to this routine has changed since earlier releases of the toolbox:
At Mark 22: |
nwkjac and njcpvt were removed from the interface; neq was made optional |
nag_ode_ivp_stiff_imp_bandjac (d02nh) is a general purpose function for integrating the initial value problem for a stiff system of implicit ordinary differential equations coupled with algebraic equations, written in the form
It is designed specifically for the case where the resulting Jacobian is a banded matrix (see the description of
jac).
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_bandjac (d02nh) 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] = d02nh(...);
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_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_bandjac (d02nh). The integrator diagnostic function
nag_ode_ivp_stiff_integ_diag (d02ny) may be called after the call to
nag_ode_ivp_stiff_imp_bandjac (d02nh). 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_bandjac (d02nh) 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_bandjac (d02nh) 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.
In general, you are advised to choose the BDF option (setup function
nag_ode_ivp_stiff_bdf (d02nv)) but if efficiency is of great importance and especially if it is suspected that
$\frac{\partial}{\partial y}\left({A}^{-1}g\right)$ has complex eigenvalues near the imaginary axis for some part of the integration, you should try the BLEND option (setup function
nag_ode_ivp_stiff_blend (d02nw)).
This example solves the well-known stiff Robertson problem written as an implicit differential system and in implicit form
exploiting the fact that we can show that
${\left(a+b+c\right)}^{\prime}=0$ for all time. Integration is over the range
$\left[0,10.0\right]$ with initial conditions
$a=1.0$ and
$b=c=0.0$ using scalar relative error control and vector absolute error control (
${\mathbf{itol}}=2$). We integrate using a BDF method (setup function
nag_ode_ivp_stiff_bdf (d02nv)) and a modified Newton method. The Jacobian is calculated numerically and we employ a default monitor, string
nag_ode_ivp_stiff_exp_fulljac_dummy_monit (d02nby). We perform a normal integration (
${\mathbf{itask}}=1$) to obtain the value at
${\mathbf{tout}}=10.0$ by integrating past
tout and interpolating. We also illustrate the use of
${\mathbf{itask}}=6$ to calculate initial values of
$y$ and
${y}^{\prime}$ and then return without integrating further.
function d02nh_example
fprintf('d02nh example results\n\n');
neq = int64(3);
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( ...
neq, sdysav, maxord, 'Newton', petzld, ...
const, tcrit, hmin, hmax, h0, maxstp, ...
mxhnil, 'Average-L2', rwork);
ml = int64(1);
mu = int64(2);
nwkjac = int64(neq*(2*ml+mu+1));
[rwork, ifail] = d02nt( ...
neq, neq, 'Numerical', ml, mu, nwkjac, neq, rwork);
t = 0;
tout = 10;
y = [1; 0; 0];
ydot = zeros(neq, 1);
rtol = [0.0001];
atol = [1e-06; 1e-07; 1e-06];
itol = int64(2);
inform(1:23) = int64(0);
ysave = zeros(neq, sdysav);
wkjac = zeros(nwkjac, 1);
jacpvt(1:neq) = int64(0);
lderiv = [false; true];
itask = int64(6);
itrace = int64(0);
[t, tout, y, ydot, rwork, inform, ysave, wkjac, jacpvt, lderiv, ifail] = ...
d02nh( ...
t, tout, y, ydot, rwork, rtol, atol, itol, inform, @resid, ysave, ...
'd02nhz', wkjac, jacpvt, 'd02nby', lderiv, itask, itrace);
fprintf ('Initial y : %8.4f %8.4f %8.4f\n',y);
fprintf ('Initial ydot: %8.4f %8.4f %8.4f\n',ydot);
itask = int64(1);
tout = 10;
[t, tout, y, ydot, rwork, inform, ysave, wkjac, jacpvt, lderiv, ifail] = ...
d02nh( ...
t, tout, y, ydot, rwork, rtol, atol, itol, inform, @resid, ysave, ...
'd02nhz', wkjac, jacpvt, 'd02nby', lderiv, itask, itrace);
fprintf ('\n x y\n');
fprintf ('%8.3f %13.5f %13.5f %13.5f\n',tout,y);
[hu, h, tcur, tolsf, nst, nre, nje, nqu, nq, niter, imxer, algequ, ifail] = ...
d02ny( ...
neq, neq, 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);
function [f, ires] = fcn(neq, t, y, ires)
f = zeros(3,1);
f(1) = -0.04*y(1) + 10000*y(2)*y(3);
f(2) = 0.04*y(1) - 10000*y(2)*y(3) - 30000000*y(2)*y(2);
f(3) = 30000000*y(2)*y(2);
function p = jac(neq, t, y, h, d, p)
p = zeros(neq, neq);
hxd = h*d;
p(1,1) = 1 - hxd*(-0.04);
p(1,2) = -hxd*(10000*y(3));
p(1,3) = -hxd*(10000*y(2));
p(2,1) = -hxd*(0.04);
p(2,2) = 1 - hxd*(-10000*y(3)-60000000*y(2));
p(2,3) = -hxd*(-10000*y(2));
p(3,2) = -hxd*(60000000*y(2));
p(3,3) = 1 - hxd*(0);
function [r, ires] = resid(neq, t, y, ydot, ires)
r = zeros(neq,1);
r(1) = -ydot(1) - ydot(2) - ydot(3);
r(2) = -ydot(2);
r(3) = -ydot(3);
if (ires == 1) ;
r(1) = 0.0d0 + r(1);
r(2) = 0.04d0*y(1) - 1.0d4*y(2)*y(3) - 3.0d7*y(2)*y(2) + r(2);
r(3) = 3.0d7*y(2)*y(2) + r(3);
end ;