Cases prefixed with W are classified as warnings and do not generate an error of type NAG:error_n. See nag_issue_warnings.
Open in the MATLAB editor: d02za_example
function d02za_example fprintf('d02za example results\n\n'); global isol tkeep ykeep % Initialize 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 = 10; petzld = false; const = zeros(6, 1); rwork = zeros(50+4*neqmax, 1); % d02nv is a setup routine to be called prior to d02nb. [const, rwork, ifail] = d02nv( ... neqmax, sdysav, maxord, 'Newton', petzld, ... const, tcrit, hmin, hmax, h0, maxstp, ... mxhnil, 'Average-L2', rwork); % d02ns determines how d02nb evaluates the Jacobian. [rwork, ifail] = d02ns( ... neq, neqmax, 'Analytical', nwkjac, rwork); % Initialize integration variables and arrays. inform(1:23) = int64(0); ysave = zeros(neq, sdysav); wkjac = zeros(nwkjac, 1); % First case. Integrate to tout without passing tout (tcrit=tout and itask=4) % use B.D.F formulae with a Newton method. % Evaluate Jacobian numerically (d02nbz), no monitoring (d02nby). t = 0.0; tout = 10.0; itask = int64(4); itrace = int64(0); y = [1; 0; 0]; itol = int64(1); rtol = [0.0001]; atol = [1e-07]; isol = 1; tkeep(isol) = t; ykeep(1:3,isol) = y; % Output initial and final solutions. fprintf(' Using Analytical Jacobian\n\n'); fprintf(' x y_1 y_2 y_3\n'); fprintf(' %8.3f %5.1f %5.1f %5.1f\n', t, y); [t, y, ydot, rwork, inform, ysave, wkjac, ifail] = ... d02nb( ... t, tout, y, rwork, rtol, atol, itol, inform, @fcn, ysave, ... @jac, wkjac, @monitr, itask, itrace); fprintf(' %8.3f %5.1f %5.1f %5.1f\n', t, y); % Plot results. fig1 = figure; display_plot(tkeep, ykeep) function [f, ires] = fcn(neq, t, y, ires) % Evaluate derivative vector. f = zeros(3,1); f(1) = -0.04*y(1) + 1.0d4*y(2)*y(3); f(2) = 0.04*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, p) % Evaluate the Jacobian. p = zeros(neq, neq); hxd = h*d; p(1,1) = 1 - hxd*(-0.04); p(1,2) = - hxd*(1.0d4*y(3)); p(1,3) = - hxd*(1.0d4*y(2)); p(2,1) = - hxd*(0.04); p(2,2) = 1 - hxd*(-1.0d4*y(3)-6.0d7*y(2)); p(2,3) = - hxd*(-1.0d4*y(2)); p(3,2) = - hxd*(6.0d7*y(2)); p(3,3) = 1 - hxd*(0); function [hnext, y, imon, inln, hmin, hmax] = ... monitr(neq, neqmax, t, hlast, hnext, y, ydot, ... ysave, r, acor, imon, hmin, hmax, nqu) global isol tkeep ykeep inln=int64(0); if (imon == int64(1)) isol = isol + 1; tkeep(isol) = t; ykeep(1:3,isol) = y; end function display_plot(x,y) % Formatting for title and axis labels. % Plot one of the curves and then add the other two. hline3 = plot(x, y(3,:)); hold on; [haxes, hline1, hline2] = plotyy(x, 100*y(2,:), x, y(1,:)); % Set the axis limits and the tick specifications to beautify the plot. set(haxes(1), 'YLim', [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.005]); 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'); % Add title. th = title({'Stiff Robertson Problem','Using BDF with Modified Newton'}); set(th,'position',[0.05,0.004]); % Label the x axis, and both y axes. xlabel('x'); ylabel(haxes(1), 'Solution (100*b,c)'); ylabel(haxes(2), 'Solution (a)'); % Add a legend. legend('c', '100*b', 'a', 'Location', 'East'); % Set some features of the three lines set(hline1, 'Linewidth', 0.5, 'Marker', '+', 'LineStyle', '-', ... 'Color', 'Magenta'); set(hline2, 'Linewidth', 0.5, 'Marker', '*', 'LineStyle', ':'); set(hline3, 'Linewidth', 0.5, 'Marker', 'x', 'LineStyle', '--'); hold off;
d02za example results
Using Analytical Jacobian
x y_1 y_2 y_3
0.000 1.0 0.0 0.0
10.000 0.8 0.0 0.2