Both interval and step oriented modes of operation are available and also modes designed to permit intermediate output within an interval oriented mode.
function d02nd_example
fprintf('d02nd example results\n\n');
global ncall ykeep tkeep tout
global xout rwork
neq = int64(3);
neqmax = int64(neq);
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);
nelts = 8;
nia = int64(neq+1);
nja = int64(nelts);
ia = int64([1; 3; 6; 9]);
ja = int64([1; 2; 1; 2; 3; 1; 2; 3]);
njcpvt = int64(150);
nwkjac = int64(50);
eta = 1.0e-4;
u = 0.1;
sens = 0.0;
lblock = true;
[jacpvt, rwork, ifail] = ...
d02nu( ...
neq, neqmax, 'Numerical', nwkjac, ia, ja, ...
njcpvt, sens, u, eta, lblock, rwork, 'nia', nia, 'nja', nja);
inform(1:23) = int64(0);
wkjac = zeros(nwkjac, 1);
ysave = zeros(neq, sdysav);
t = 0;
tout = 10.0;
itask = int64(1);
itrace = int64(0);
y = [1; 0; 0];
itol = int64(1);
rtol = [1e-4];
atol = [1e-7];
isplit = int64(0);
fprintf('Numerical Jacobian, structure not supplied\n\n');
fprintf(' x y_1 y_2 y_3\n');
fprintf('%8.3f %8.5f %8.5f %8.5f\n',t,y);
xout = 2.0;
[t, y, ydot, rwork, inform, ysave, wkjac, jacpvt,ifail] = ...
d02nd( ...
t, tout, y, rwork, rtol, atol, itol, inform, @fcn, ysave, ...
'd02ndz', wkjac, jacpvt, @monitr1, itask, itrace);
[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);
t = 0;
y = [1; 0; 0];
isplit = int64(0);
ncall = int64(1);
tkeep = t;
ykeep = y;
[const, rwork, ifail] = d02nv( ...
neqmax, sdysav, maxord, 'Newton', petzld, ...
const, tcrit, hmin, hmax, h0, maxstp, ...
mxhnil, 'Average-L2', rwork);
[jacpvt, rwork, ifail] = d02nu( ...
neq, neqmax, 'Structural', nwkjac, ia, ja, ...
njcpvt, sens, u, eta, lblock, rwork, 'nia', ...
nia, 'nja', nja);
fprintf('Numerical Jacobian, structure supplied\n\n');
fprintf(' x y_1 y_2 y_3\n');
fprintf('%8.3f %8.5f %8.5f %8.5f\n',t,y);
xout = 2.0;
[t, y, ydot, rwork, inform, ysave, wkjac, jacpvt,ifail] = ...
d02nd( ...
t, tout, y, rwork, rtol, atol, itol, inform, @fcn, ysave, ...
'd02ndz', wkjac, jacpvt, @monitr2, itask, itrace);
[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);
[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);
fig1 = figure;
display_plot(tkeep,ykeep)
function [hnext, y, imon, inln, hmin, hmax] = ...
monitr1(neq, neqmax, t, hlast, hnext, y, ydot, ysave, ...
r, acor, imon, hmin, hmax, nqu)
[imon] = interp(neq, t, hlast, hnext, ysave, imon, nqu);
inln = int64(0);
function [hnext, y, imon, inln, hmin, hmax] = ...
monitr2( neq, neqmax, t, hlast, hnext, y, ydot, ysave, ...
r, acor, imon, hmin, hmax, nqu)
global ncall ykeep tkeep tout
if (imon == 1 && t>0.00001 && t<tout)
ncall = ncall + 1;
ykeep(:,ncall) = y;
tkeep (ncall) = t;
end
[imon] = interp(neq, t, hlast, hnext, ysave, imon, nqu);
inln = int64(0);
function [imon] = interp(neq, t, hlast, hnext, ysave, imon, nqu)
global xout rwork
dims = size(ysave);
lbeg = dims(1) + 50 + 1;
lend = lbeg + neq - 1;
lw(1:neq) = rwork(lbeg:lend);
while (xout <= 10.0 && t-hlast < xout && xout <= t)
[sol, ifail] = d02xk(xout, neq, ysave, lw, t, nqu, hlast, hnext);
fprintf('%8.3f %13.5f %13.5f %13.5f\n', xout, sol(1:3));
xout = xout + 2;
end
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 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 1.2]);
set(haxes(1), 'XMinorTick', 'on', 'YMinorTick', 'on');
set(haxes(1), 'YTick', [0:0.2:1]);
set(haxes(2), 'YLim', [0 4e-5]);
set(haxes(2), 'YMinorTick', 'on');
set(haxes(2), 'YTick', [5e-6:5e-6:3.5e-5]);
set(haxes(1), 'XLim', [-0.1 10]);
set(haxes(2), 'XLim', [-0.1 10]);
ht = title({'Stiff Robertson Problem: BDF Method', ...
'Jacobian: Sparse with Structure Supplied'});
set(ht,'Position',[5,1.05,0.0]);
xlabel('x');
ylabel(haxes(1),'Solution (a,c)');
ylabel(haxes(2),'Solution (b)');
legend('c','a','b','Location','West');
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;