y^{′} = f(t,y) given y(t_{0}) = y_{0}
$${y}^{\prime}=f(t,y)\text{\hspace{1em} given \hspace{1em}}y\left({t}_{0}\right)={y}_{0}$$ |
None.
None.
Cases prefixed with W are classified as warnings and do not generate an error of type NAG:error_n. See nag_issue_warnings.
y^{ ′ ′ } = − y, y(0) = 0, y^{′}(0) = 1
$${y}^{\prime \prime}=-y\text{, \hspace{1em}}y\left(0\right)=0\text{, \hspace{1em}}{y}^{\prime}\left(0\right)=1$$ |
y_{1}^{ ′ } = y_{2}
$${y}_{1}^{\prime}={y}_{2}$$ |
y_{2}^{ ′ } = − y_{1}
$${y}_{2}^{\prime}=-{y}_{1}$$ |
Open in the MATLAB editor: nag_ode_ivp_rk_range_example
function nag_ode_ivp_rk_range_example % Initialize variables and arrays. tstart = 0; ystart = [0; 1]; hstart = 0; tend = 2*pi; tol = 0.001; thres = [1e-08; 1e-08]; method = 1; task = 'Usual Task'; errass = false; lenwrk = 64; neq = 2; twant = 0.7853981633974483; ygot = [0; 0]; ymax = [0; 0]; fprintf('nag_ode_ivp_rk_range example program results \n\n'); % We run through the calculation twice; the first time with a small % number of points (when we'll output the results), and the second time % with a larger number (when we'll plot them). npts = [8, 40]; for ipts = 1:2 tinc = (tend-tstart)/npts(ipts); ncall = [0, 0]; y1keep = [ystart(1),ystart(1)]; y2keep = [ystart(2),ystart(2)]; err = [0, 0]; tolkeep = [1]; tkeep = [tstart]; for i = 1:2 if (i == 1) tol = double(1e-3); elseif (i == 2) tol = double(1e-4); end tolkeep(i) = tol; % nag_ode_ivp_rk_setup is a setup routine to be called prior to nag_ode_ivp_rk_range. [work, ifail] = nag_ode_ivp_rk_setup(tstart, ystart, tend, tol, thres, ... int64(method), task, errass, int64(lenwrk), 'neq', ... int64(neq), 'hstart', hstart); if ifail ~= 0 % Illegal input. Print message and exit. error('Warning: nag_ode_ivp_rk_setup returned with ifail = %1d ',ifail); end if ipts == 1 fprintf('Calculation with tol = %1.3e\n\n',tol); disp(' t y1 y2'); fprintf(' %1.3f %9.3f %9.3f\n', tstart, ystart(1), ystart(2)); end for j = npts(ipts)-1:-1:0 twant = tend - j*tinc; [tgot, ygot, ypgot, ymax, work, ifail] = ... nag_ode_ivp_rk_range(@f, int64(neq), twant, ygot, ymax, work); if ifail ~= 0 % Illegal input. Print message and exit. error('Warning: nag_ode_ivp_rk_range returned with ifail = %1d ',... ifail); end if ipts == 1 fprintf(' %1.3f %9.3f %9.3f\n', tgot, ygot(1), ygot(2)); end % Store results for plotting (second time around). ncall(i) = ncall(i) + 1; y1keep(ncall(i),i) = ygot(1); y2keep(ncall(i),i) = ygot(2); err(ncall(i),i) = ygot(1)-sin(tgot); tkeep(ncall(i),i) = tgot; end % nag_ode_ivp_rk_diag is a diagnostic routine. [totfcn, stpcst, waste, stpsok, hnext, ifail] = nag_ode_ivp_rk_diag; if ipts == 1 fprintf(['\nCost of the integration in evaluations of ', ... 'F is %1.0f\n\n\n'], totfcn); end end end % Plot results. fig = figure('Number', 'off'); display_plot(tkeep, y1keep, y2keep, err); function [yp] = f(t, y) % Evaluate derivative vector. yp = zeros(2, 1); yp(1) = y(2); yp(2) = -y(1); function display_plot(tkeep, y1keep, y2keep, error1) % Formatting for title and axis labels. titleFmt = {'FontName', 'Helvetica', 'FontWeight', 'Bold', 'FontSize', 14}; labFmt = {'FontName', 'Helvetica', 'FontWeight', 'Bold', 'FontSize', 13}; set(gca, 'FontSize', 13); % for legend, axis tick labels, etc. % Plot two of the curves and then add the other one. [haxes, hline1, hline2] = plotyy(tkeep(:,1), y1keep(:,2), tkeep(:,1), ... abs(error1(:,1)), 'plot', 'semilogy'); hold on hline3 = plot(haxes(1), tkeep(:,1), y2keep(:,2)); % Select the required y axis, and plot the fourth curve. axes(haxes(2)); hold on hline4 = plot(tkeep(:,1), abs(error1(:,2))); % Set the axis limits and the tick specifications to beautify the plot. set(haxes(1), 'YLim', [-1 1.2]); set(haxes(1), 'YMinorTick', 'on'); set(haxes(1), 'YTick', [-1 -0.6 -0.2 0.2 0.6 1.0]); set(haxes(2), 'YLim', [1e-7 1e-2]); set(haxes(2), 'YMinorTick', 'on'); set(haxes(2), 'YTick', [1e-7 1e-6 1e-5 1e-4 1e-3 1e-2]); for iaxis = 1:2 % These properties must be the same for both sets of axes. set(haxes(iaxis), 'XLim', [0 7]); set(haxes(iaxis), 'XTick', [0 1 2 3 4 5 6 7]); set(haxes(iaxis), 'FontSize', 13); end set(gca, 'box', 'off'); % so ticks aren't shown on opposite axes. % Add title. title(['1st-order ODEs using RK ', ... 'Low-order Method with 2 Tolerances'], titleFmt{:}) % Label the x axis, and both y axes. xlabel('t', labFmt{:}); ylabel(haxes(1),'Solution (y,y'')', labFmt{:}); ylabel(haxes(2),'abs(Error)', labFmt{:}); % Add a legend. legend('y-error (tol= 0.001)','y-error (tol= 0.0001)','y','y''', ... 'Location','Best') % Set some features of the lines. set(hline1, 'Linewidth', 0.5, 'Marker', '+','LineStyle','-'); set(hline2, 'Linewidth', 0.5, 'Marker', 's','LineStyle','--'); set(hline3, 'Linewidth', 0.5, 'Marker', 'x','LineStyle','-.'); set(hline4, 'Linewidth', 0.5, 'Marker', '*','LineStyle',':');
nag_ode_ivp_rk_range example program results Calculation with tol = 1.000e-03 t y1 y2 0.000 0.000 1.000 0.785 0.707 0.707 1.571 0.999 -0.000 2.356 0.706 -0.706 3.142 -0.000 -0.999 3.927 -0.706 -0.706 4.712 -0.998 0.000 5.498 -0.705 0.706 6.283 0.001 0.997 Cost of the integration in evaluations of F is 124 Calculation with tol = 1.000e-04 t y1 y2 0.000 0.000 1.000 0.785 0.707 0.707 1.571 1.000 -0.000 2.356 0.707 -0.707 3.142 -0.000 -1.000 3.927 -0.707 -0.707 4.712 -1.000 0.000 5.498 -0.707 0.707 6.283 0.000 1.000 Cost of the integration in evaluations of F is 235
Open in the MATLAB editor: d02pc_example
function d02pc_example % Initialize variables and arrays. tstart = 0; ystart = [0; 1]; hstart = 0; tend = 2*pi; tol = 0.001; thres = [1e-08; 1e-08]; method = 1; task = 'Usual Task'; errass = false; lenwrk = 64; neq = 2; twant = 0.7853981633974483; ygot = [0; 0]; ymax = [0; 0]; fprintf('d02pc example program results \n\n'); % We run through the calculation twice; the first time with a small % number of points (when we'll output the results), and the second time % with a larger number (when we'll plot them). npts = [8, 40]; for ipts = 1:2 tinc = (tend-tstart)/npts(ipts); ncall = [0, 0]; y1keep = [ystart(1),ystart(1)]; y2keep = [ystart(2),ystart(2)]; err = [0, 0]; tolkeep = [1]; tkeep = [tstart]; for i = 1:2 if (i == 1) tol = double(1e-3); elseif (i == 2) tol = double(1e-4); end tolkeep(i) = tol; % d02pv is a setup routine to be called prior to d02pc. [work, ifail] = d02pv(tstart, ystart, tend, tol, thres, ... int64(method), task, errass, int64(lenwrk), 'neq', ... int64(neq), 'hstart', hstart); if ifail ~= 0 % Illegal input. Print message and exit. error('Warning: d02pv returned with ifail = %1d ',ifail); end if ipts == 1 fprintf('Calculation with tol = %1.3e\n\n',tol); disp(' t y1 y2'); fprintf(' %1.3f %9.3f %9.3f\n', tstart, ystart(1), ystart(2)); end for j = npts(ipts)-1:-1:0 twant = tend - j*tinc; [tgot, ygot, ypgot, ymax, work, ifail] = ... d02pc(@f, int64(neq), twant, ygot, ymax, work); if ifail ~= 0 % Illegal input. Print message and exit. error('Warning: d02pc returned with ifail = %1d ',... ifail); end if ipts == 1 fprintf(' %1.3f %9.3f %9.3f\n', tgot, ygot(1), ygot(2)); end % Store results for plotting (second time around). ncall(i) = ncall(i) + 1; y1keep(ncall(i),i) = ygot(1); y2keep(ncall(i),i) = ygot(2); err(ncall(i),i) = ygot(1)-sin(tgot); tkeep(ncall(i),i) = tgot; end % d02py is a diagnostic routine. [totfcn, stpcst, waste, stpsok, hnext, ifail] = d02py; if ipts == 1 fprintf(['\nCost of the integration in evaluations of ', ... 'F is %1.0f\n\n\n'], totfcn); end end end % Plot results. fig = figure('Number', 'off'); display_plot(tkeep, y1keep, y2keep, err); function [yp] = f(t, y) % Evaluate derivative vector. yp = zeros(2, 1); yp(1) = y(2); yp(2) = -y(1); function display_plot(tkeep, y1keep, y2keep, error1) % Formatting for title and axis labels. titleFmt = {'FontName', 'Helvetica', 'FontWeight', 'Bold', 'FontSize', 14}; labFmt = {'FontName', 'Helvetica', 'FontWeight', 'Bold', 'FontSize', 13}; set(gca, 'FontSize', 13); % for legend, axis tick labels, etc. % Plot two of the curves and then add the other one. [haxes, hline1, hline2] = plotyy(tkeep(:,1), y1keep(:,2), tkeep(:,1), ... abs(error1(:,1)), 'plot', 'semilogy'); hold on hline3 = plot(haxes(1), tkeep(:,1), y2keep(:,2)); % Select the required y axis, and plot the fourth curve. axes(haxes(2)); hold on hline4 = plot(tkeep(:,1), abs(error1(:,2))); % Set the axis limits and the tick specifications to beautify the plot. set(haxes(1), 'YLim', [-1 1.2]); set(haxes(1), 'YMinorTick', 'on'); set(haxes(1), 'YTick', [-1 -0.6 -0.2 0.2 0.6 1.0]); set(haxes(2), 'YLim', [1e-7 1e-2]); set(haxes(2), 'YMinorTick', 'on'); set(haxes(2), 'YTick', [1e-7 1e-6 1e-5 1e-4 1e-3 1e-2]); for iaxis = 1:2 % These properties must be the same for both sets of axes. set(haxes(iaxis), 'XLim', [0 7]); set(haxes(iaxis), 'XTick', [0 1 2 3 4 5 6 7]); set(haxes(iaxis), 'FontSize', 13); end set(gca, 'box', 'off'); % so ticks aren't shown on opposite axes. % Add title. title(['1st-order ODEs using RK ', ... 'Low-order Method with 2 Tolerances'], titleFmt{:}) % Label the x axis, and both y axes. xlabel('t', labFmt{:}); ylabel(haxes(1),'Solution (y,y'')', labFmt{:}); ylabel(haxes(2),'abs(Error)', labFmt{:}); % Add a legend. legend('y-error (tol= 0.001)','y-error (tol= 0.0001)','y','y''', ... 'Location','Best') % Set some features of the lines. set(hline1, 'Linewidth', 0.5, 'Marker', '+','LineStyle','-'); set(hline2, 'Linewidth', 0.5, 'Marker', 's','LineStyle','--'); set(hline3, 'Linewidth', 0.5, 'Marker', 'x','LineStyle','-.'); set(hline4, 'Linewidth', 0.5, 'Marker', '*','LineStyle',':');
d02pc example program results Calculation with tol = 1.000e-03 t y1 y2 0.000 0.000 1.000 0.785 0.707 0.707 1.571 0.999 -0.000 2.356 0.706 -0.706 3.142 -0.000 -0.999 3.927 -0.706 -0.706 4.712 -0.998 0.000 5.498 -0.705 0.706 6.283 0.001 0.997 Cost of the integration in evaluations of F is 124 Calculation with tol = 1.000e-04 t y1 y2 0.000 0.000 1.000 0.785 0.707 0.707 1.571 1.000 -0.000 2.356 0.707 -0.707 3.142 -0.000 -1.000 3.927 -0.707 -0.707 4.712 -1.000 0.000 5.498 -0.707 0.707 6.283 0.000 1.000 Cost of the integration in evaluations of F is 235