At Mark 22: | njcpvt was removed from the interface; neq was made optional |
$$A\left(t,y\right){y}^{\prime}=g\left(t,y\right)\text{.}$$ |
% % data initialization % call linear algebra setup routine call integrator setup routine irevcm = int32(0); [..., irevcm, ...] = d02nn(); while (irevcm > 0) if (irevcm > 7 and irevcm < 11) if (irevcm == 8) supply the jacobian matrix (i) elseif (irevcm == 9) perform monitoring tasks requested by the user (ii) elseif (irevcm == 10) indicates an unsuccessful step end else evaluate the residual (iii) end [..., irevcm, ...] = d02nn(); end % % post processing (optional linear algebra diagnostic call % (sparse case only), optional integrator diagnostic call) %
(i) | Supply the Jacobian matrix
You need only provide this facility if the argument ${\mathbf{jceval}}=\text{'A'}$ (or ${\mathbf{jceval}}=\text{'F'}$ if using sparse matrix linear algebra) in a call to the linear algebra setup function (see jceval in nag_ode_ivp_stiff_sparjac_setup (d02nu)). If the Jacobian matrix is to be evaluated numerically by the integrator, then the remainder of section (i) can be ignored. We must define the system of nonlinear equations which is solved internally by the integrator. The time derivative, ${y}^{\prime}$, has the form
The system of nonlinear equations that is solved has the form
It is the Jacobian matrix $\frac{\partial r}{\partial y}$ that you must supply as follows:
Hereafter in this document this operation will be referred to as jac. |
|||||||||||||||||||||||
(ii) | Perform tasks requested by you
This operation is essentially a monitoring function and additionally provides the opportunity of changing the current values of y, ydot, hnext (the step size that the integrator proposes to take on the next step), hmin (the minimum step size to be taken on the next step), and hmax (the maximum step size to be taken on the next step).
The scaled local error at the end of a time step may be obtained by calling the double function nag_ode_ivp_stiff_errest (d02za) as follows:
[errloc, ifail] = d02za(rwork(51+neqmax:51+neqmax+neq-1), rwork(51:51+neq-1)); % Check ifail before proceeding The following gives details of the location within the array rwork of variables that may be of interest to you:
You are advised to consult the description of
monitr in nag_ode_ivp_stiff_imp_fulljac (d02ng)
for details on what optional input can be made.
If either y or ydot are changed, then imon must be set to $2$ before return to nag_ode_ivp_stiff_imp_revcom (d02nn). If either of the values hmin or hmax are changed, then imon must be set $\text{}\ge 3$ before return to nag_ode_ivp_stiff_imp_revcom (d02nn). If hnext is changed, then imon must be set to $4$ before return to nag_ode_ivp_stiff_imp_revcom (d02nn).
In addition you can force nag_ode_ivp_stiff_imp_revcom (d02nn) to evaluate the residual vector
Hereafter in this document this operation will be referred to as monitr. |
|||||||||||||||||||||||
(iii) | Evaluate the residual
This operation must evaluate the residual
Hereafter in this document this operation will be referred to as resid. |
itol | rtol | atol | ${w}_{i}$ |
1 | scalar | scalar | ${\mathbf{rtol}}\left(1\right)\times \left|{y}_{i}\right|+{\mathbf{atol}}\left(1\right)$ |
2 | scalar | vector | ${\mathbf{rtol}}\left(1\right)\times \left|{y}_{i}\right|+{\mathbf{atol}}\left(i\right)$ |
3 | vector | scalar | ${\mathbf{rtol}}\left(i\right)\times \left|{y}_{i}\right|+{\mathbf{atol}}\left(1\right)$ |
4 | vector | vector | ${\mathbf{rtol}}\left(i\right)\times \left|{y}_{i}\right|+{\mathbf{atol}}\left(i\right)$ |
[j, iplace] = d02nr(inform);
Cases prefixed with W are classified as warnings and do not generate an error of type NAG:error_n. See nag_issue_warnings.
$$\begin{array}{l}{r}_{1}=\left({a}^{\prime}+{b}^{\prime}+{c}^{\prime}\right)\\ {r}_{2}=0.04a-1.0\mathrm{E}4bc-3.0\mathrm{E}7{b}^{2}-{b}^{\prime}\\ {r}_{3}=\phantom{0.04a-1.0\mathrm{E}4bc-}3.0\mathrm{E}7{b}^{2}-{c}^{\prime}\end{array}$$ |
Open in the MATLAB editor: d02nn_example
function d02nn_example fprintf('d02nn 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 = 0; 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); % Sparse Jacobian supplied, setup ia = int64(0); ja = int64(0); njcpvt = int64(150); nwkjac = int64(100); eta = 1.0e-4; u = 0.1; sens = 0.0; lblock = true; [jacpvt, rwork, ifail] = d02nu(... neq, neqmax, 'Analytical', nwkjac, ia, ja, ... njcpvt, sens, u, eta, lblock, rwork); % Integration input variable intialization t = 0; tout = 10; y = [1; 0; 0]; ydot = [0; 0; 0]; rtol = [0.0001]; atol = [1e-07]; itol = int64(1); inform = zeros(23, 1, 'int64'); ysave = zeros(neq, sdysav); wkjac = zeros(nwkjac,1); imon = int64(0); inln = int64(0); ires = int64(1); irevcm = int64(0); lderiv = [false; false]; itask = int64(3); itrace = int64(0); nfails = 0; % pointers into rwork locations lacorb = neqmax + 50; lsavrb = lacorb + neqmax; l1 = lsavrb+1; l2 = l1+1; l3 = l2+1; m1 = lacorb+1; m2 = m1+1; m3 = m2+1; fprintf(' Analytic Jacobian\n\n'); fprintf(' x y_1 y_2 y_3\n'); fprintf(' %8.3f %5.1f %5.1f %5.1f\n', t, y); first_time = true; % Main reverse communication loop controlled by irevcm while (irevcm ~= 0 || first_time) first_time = false; [t, tout, y, ydot, rwork, inform, ysave, wkjac, ... jacpvt, imon, inln, ires, irevcm, lderiv, ifail] = ... d02nn(t, tout, y, ydot, rwork, rtol, atol, itol, inform, ... ysave, wkjac, jacpvt, imon, inln, ires, irevcm, lderiv, itask, itrace); if irevcm == 1 || irevcm == 3 || irevcm == 6 || irevcm == 11 % Equivalent to resid evaluation in forward communication routines % ydot stored in ydot, resid returned in rwork(l1) rwork(l1) = -ydot(1) - ydot(2) - ydot(3); rwork(l2) = -ydot(2); rwork(l3) = -ydot(3); if (ires == 1) rwork(l1) = rwork(l1); rwork(l2) = 0.04*y(1) - 1e4*y(2)*y(3) - 3e7*y(2)*y(2) + rwork(l2); rwork(l3) = 3e7*y(2)*y(2) + rwork(l3); end elseif (irevcm == 2) % Equivalent to resid evaluation in forward communication routines % ydot stored in rwork(l1:), resid returned in rwork(l1) rwork(l1) = -rwork(l1) - rwork(l2) - rwork(l3); rwork(l2) = -rwork(l2); rwork(l3) = -rwork(l3); elseif (irevcm == 4 || irevcm == 7) % Equivalent to resid evaluation in forward communication routines % ydot stored in ydot, resid returned in rwork(m1) rwork(m1) = -ydot(1) - ydot(2) - ydot(3); rwork(m2) = -ydot(2); rwork(m3) = -ydot(3); if (ires == 1) rwork(m1) = rwork(m1); rwork(m2) = 0.04*y(1) - 1e4*y(2)*y(3) - 3e7*y(2)*y(2) + rwork(m2); rwork(m3) = 3e7*y(2)*y(2) + rwork(m3); end elseif (irevcm == 5) % Equivalent to resid evaluation in forward communication routines % ydot stored in rwork(l1:), resid returned in ydot ydot(1) = -rwork(l1) - rwork(l2) - rwork(l); ydot(2) = 0.04*y(1) - 1e4*y(2)*y(3) - 3e7*y(2)*y(2) - rwork(l2); ydot(3) = 3e7*y(2)*y(2) - rwork(l3); elseif (irevcm == 8) % Equivalent to jac evaluation in forward communication routines [j, iplace] = d02nr(inform); hxd = rwork(16)*rwork(20); if (iplace < 2) if (j < 2) rwork(l1) = 1 - hxd*(0); rwork(l2) = 0 - hxd*(0.04); rwork(l3) = 0 - hxd*(0); elseif (j == 2) rwork(l1) = 1 - hxd*(0); rwork(l2) = 1 - hxd*(-1e4*y(3)-6e7*y(2)); rwork(l3) = 0 - hxd*(6e7*y(2)); elseif (j > 2) rwork(l1) = 1 - hxd*(0); rwork(l2) = 0 - hxd*(-1e4*y(2)); rwork(l3) = 1 - hxd*(0); end else if (j < 2) rwork(m1) = 1 - hxd*(0); rwork(m2) = 0 - hxd*(0.04); rwork(m3) = 0 - hxd*(0); elseif (j == 2) rwork(m1) = 1 - hxd*(0); rwork(m2) = 1 - hxd*(-1e4*y(3)-6e7*y(2)); rwork(m3) = 0 - hxd*(6e7*y(2)); elseif (j > 2) rwork(m1) = 1 - hxd*(0); rwork(m2) = 0 - hxd*(-1e4*y(2)); rwork(m3) = 1 - hxd*(0); end end elseif (irevcm == 10) % Step failure nfails = nfails + 1; end end [hu, h, tcur, tolsf, nst, nre, nje, nqu, nq, nit, imxer, algequ, ifail] = ... d02ny(neq, neqmax, rwork, inform); [y, ifail] = d02mz(tout, neq, neq, ysave, rwork); % Print solution and diagnostics fprintf(' %8.3f %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); fprintf(' number of failed steps = %5d\n\n',nfails); icall = int64(0); [liwreq, liwusd, lrwreq, lrwusd, nlu, nz, ngp, isplit, igrow, nblock] = ... d02nx(icall, true, 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);
d02nn example results Analytic Jacobian x y_1 y_2 y_3 0.000 1.0 0.0 0.0 Warning: Equation(=i1) and possibly other equations are implicit and in calculating the initial values the equations will be treated as implicit. In above message i1 = 1 10.767 0.8 0.0 0.2 Diagnostic information integration status: last and next step sizes = 0.90181, 0.90181 integration stopped at x = 10.76669 algorithm statistics: number of time-steps and Newton iterations = 55 80 number of residual and jacobian evaluations = 89 17 order of method last used and next to use = 4 4 component with largest error = 3 number of failed steps = 4 sparse jacobian statistics: njcpvt - required = 99, used = 150 nwkjac - required = 31, used = 75 No. of LU-decomps = 17, No. of nonzeros = 8 No. of function calls to form Jacobian = 0, try isplit = 73 Growth estimate = 1531, No. of blocks on diagonal 1