None.
Open in the MATLAB editor: d02nz_example
function d02nz_example fprintf('d02nz example results\n\n'); % Initialize variables and arrays for setup routine n = int64(3); ord = int64(11); sdy = int64(ord + 3); petzld = false; con = zeros(6); tcrit = 0; hmin = 1.0e-10; hmax = 10; h0 = 0; maxstp = int64(200); mxhnil = int64(5); lrwork = 50 + 4*n; rwork = zeros(lrwork); [const, rwork, ifail] = d02nw(n, sdy, ord, con, tcrit, hmin, hmax, h0, ... maxstp, mxhnil, 'Average-L2', rwork); % Setup for banded Jacbian using d02nt ml = int64(1); mu = int64(2); nwkjac = int64(15); [rwork, ifail] = d02nt(n, n, 'Analytical', ml, mu, nwkjac, n, rwork); % Initialize variables and arrays for integration t = 0; tout = 5; y = [1; 0; 0]; rtol = [0.0001]; atol = [1e-07; 1e-08; 1e-07]; itol = int64(2); inform = zeros(23, 1, 'int64'); ysave = zeros(n, sdy); wkjac = zeros(nwkjac, 1); jacpvt = zeros(n, 1, 'int64'); itask = int64(1); itrace = int64(0); % Integrate ODE from t=0 to t=tout, no monitoring, using d02nc [t, y, ydot, rwork, inform, ysave, wkjac, jacpvt, ifail] = ... d02nc(t, tout, y, rwork, rtol, atol, itol, inform, @fcn, ysave, @jac, ... wkjac, jacpvt, 'd02nby', itask, itrace); fprintf('Solution y and derivative y'' at t = %7.4f is:\n',t); fprintf('\n %10s %10s\n','y','y'''); for i=1:n fprintf(' %10.4f %10.4f\n',y(i),ydot(i)); end % Call to d02nz to alow continuation of integration up to t=10. [rwork, ifail] = d02nz(n, tcrit, h0, hmin, hmax, maxstp, mxhnil, rwork); tout = 10; % Integrate ODE from t to t=tout, no monitoring, using d02nc [t, y, ydot, rwork, inform, ysave, wkjac, jacpvt, ifail] = ... d02nc(t, tout, y, rwork, rtol, atol, itol, inform, @fcn, ysave, @jac, ... wkjac, jacpvt, 'd02nby', itask, itrace); fprintf('Solution y and derivative y'' at t = %7.4f is:\n',t); fprintf('\n %10s %10s\n','y','y'''); for i=1:n fprintf(' %10.4f %10.4f\n',y(i),ydot(i)); end function [f, ires] = fcn(neq, t, y, ires) % Evaluate derivative vector. 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 p = jac(neq, t, y, h, d, ml, mu, pIn) % Evaluate the Jacobian. p = zeros(ml+mu+1, neq); hxd = h*d; p(1,1) = 1.0d0 - hxd*(-0.04d0); p(2,1) = -hxd*(1.0d4*y(3)); p(3,1) = -hxd*(1.0d4*y(2)); p(1,2) = -hxd*(0.04d0); p(2,2) = 1.0d0 - hxd*(-1.0d4*y(3)-6.0d7*y(2)); p(3,2) = -hxd*(-1.0d4*y(2)); p(1,3) = -hxd*(6.0d7*y(2)); p(2,3) = 1.0d0 - hxd*(0.0d0);
d02nz example results Solution y and derivative y' at t = 5.0000 is: y y' 0.8915 -0.0124 0.0000 -0.0000 0.1085 0.0124 Solution y and derivative y' at t = 10.0000 is: y y' 0.8414 -0.0078 0.0000 -0.0000 0.1586 0.0078