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}$$ |
(i) | The more efficient way is to step past the point where a solution is desired, and then call nag_ode_ivp_rk_interp (d02px) to get an answer there. Within the span of the current step, you can get all the answers you want at very little cost by repeated calls to nag_ode_ivp_rk_interp (d02px). This is very valuable when you want to find where something happens, e.g., where a particular solution component vanishes. You cannot proceed in this way with method = 3${\mathbf{method}}=3$. |
(ii) | The other way to get an answer at a specific point is to set tend to this value and integrate to tend. nag_ode_ivp_rk_onestep (d02pd) will not step past tend, so when a step would carry it past, it will reduce the step size so as to produce an answer at tend exactly. After getting an answer there (tnow = tend${\mathbf{tnow}}={\mathbf{tend}}$), you can reset tend to the next point where you want an answer, and repeat. tend could be reset by a call to nag_ode_ivp_rk_setup (d02pv), but you should not do this. You should use nag_ode_ivp_rk_reset_tend (d02pw) instead because it is both easier to use and much more efficient. This way of getting answers at specific points can be used with any of the available methods, but it is the only way with method = 3${\mathbf{method}}=3$. It can be inefficient. Should this be the case, the code will bring the matter to your attention. |
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_onestep_example
function nag_ode_ivp_rk_onestep_example % Initialize variables and arrays. tstart = 0; hstart = 0; ystart = [0; 1]; tend = 2*pi; thres = [1e-08; 1e-08]; method = 2; task = 'Complex Task'; errass = false; neq = 2; lenwrk = 32*neq; y1keep(1,:) = [ystart]; y2keep(1,:) = zeros(1,2); y3keep(1,:) = zeros(1,2); tkeep(1,:) = zeros(1,2); fprintf('nag_ode_ivp_rk_onestep example program results\n\n'); % We run through the calculation twice with two tolerance values; we % output and plot the results for both of them. tol = [1e-04; 1e-05]; for itol = 1:2 % nag_ode_ivp_rk_setup is a setup routine to be called prior to nag_ode_ivp_rk_onestep. [work, ifail] = nag_ode_ivp_rk_setup(tstart, ystart, tend, tol(itol), thres, ... int64(method), task, errass, int64(lenwrk), 'neq', ... int64(neq)); if ifail ~= 0 % Illegal input. Print message and exit. error('Warning: nag_ode_ivp_rk_setup returned with ifail = %1d ',ifail); end fprintf('Calculation with tol = %1.3e\n\n',tol(itol)); disp(' t y1 y1'''); fprintf('%1.3f %8.4f %8.4f\n', tstart, ystart(1), ystart(2)); % Initialize the storage counter, and the current t value before % looping over t. ncall = 0; tnow = tstart; while (tnow < tend) [tnow, ynow, ypnow, work, ifail] = nag_ode_ivp_rk_onestep(@f, ... int64(neq), work); if ifail ~= 0 % Illegal input. Print message and exit. error('Warning: nag_ode_ivp_rk_onestep returned with ifail = %1d \n\n',... ifail); end % Output results, then store them for plotting. fprintf('%1.3f %8.4f %8.4f\n', tnow, ynow(1), ynow(2)); ncall = ncall + 1; tkeep(ncall,itol) = tnow; if (itol == 2) y1keep(ncall,:) = [ynow(1) ynow(2)]; y2keep(ncall,:) = [ynow(1)-sin(tnow) ynow(2)-cos(tnow)]; else y3keep(ncall,:) = [ynow(1)-sin(tnow) ynow(2)-cos(tnow)]; end end % nag_ode_ivp_rk_diag is a diagnostic routine. [totfcn, stpcst, waste, stpsok, hnext, ifail] = nag_ode_ivp_rk_diag; fprintf(['\nCost of the integration in evaluations of F is ' ... '%1.0f\n\n'], totfcn); end % Plot results. fig = figure('Number', 'off'); display_plot(tkeep, y1keep, y2keep, y3keep) 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,y3keep) % 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. % Calculate spline through data values. xx = tkeep(1,2):0.01:tkeep(length(tkeep(:,2)),2); yy1 = spline(tkeep(:,2),y1keep(:,1),xx); yy2 = spline(tkeep(:,2),y1keep(:,2),xx); % Plot two of the curves and then add the other three. [haxes, hline1, hline2] = plotyy(xx, yy2, tkeep(1:length(y3keep),1), ... abs(y3keep(:,1)), 'plot', 'semilogy'); hold on plot(xx,yy1,':'); hold on plot(tkeep(:,2),y1keep(:,1),'+'); hold on plot(tkeep(:,2),y1keep(:,2),'x'); % Select the required y axis, and plot the sixth curve. axes(haxes(2)); hold on; hline3 = plot(tkeep(:,2), abs(y2keep(:,1))); % 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 1.2]); set(haxes(2), 'YLim', [1e-9 1e-4]); set(haxes(2), 'YMinorTick', 'on'); set(haxes(2), 'YTick', [1e-9 1e-8 1e-7 1e-6 1e-5 1e-4]); 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 Step-by-step RK ',... 'Medium-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.0001)','y-error (tol= 0.00001)','cos(x)',... 'sin(x)','y-solution','y''-solution','Location','Best'); % Set some features of the lines. set(hline1, 'Linewidth', 0.5,'LineStyle','--'); set(hline2, 'Linewidth', 0.5, 'Marker', 's','LineStyle','-.'); set(hline3, 'Linewidth', 0.5, 'Marker', '*','LineStyle',':'); set(hline3,'Visible','off');
nag_ode_ivp_rk_onestep example program results Calculation with tol = 1.000e-04 t y1 y1' 0.000 0.0000 1.0000 0.785 0.7071 0.7071 1.519 0.9987 0.0513 2.282 0.7573 -0.6531 2.911 0.2285 -0.9735 3.706 -0.5348 -0.8450 4.364 -0.9399 -0.3414 5.320 -0.8209 0.5710 5.802 -0.4631 0.8863 6.283 0.0000 1.0000 Cost of the integration in evaluations of F is 78 Calculation with tol = 1.000e-05 t y1 y1' 0.000 0.0000 1.0000 0.393 0.3827 0.9239 0.785 0.7071 0.7071 1.416 0.9881 0.1538 1.870 0.9557 -0.2943 2.204 0.8062 -0.5916 2.761 0.3711 -0.9286 3.230 -0.0880 -0.9961 3.587 -0.4304 -0.9026 4.022 -0.7710 -0.6368 4.641 -0.9974 -0.0717 5.152 -0.9049 0.4256 5.521 -0.6903 0.7235 5.902 -0.3718 0.9283 6.283 0.0000 1.0000 Cost of the integration in evaluations of F is 118
Open in the MATLAB editor: d02pd_example
function d02pd_example % Initialize variables and arrays. tstart = 0; hstart = 0; ystart = [0; 1]; tend = 2*pi; thres = [1e-08; 1e-08]; method = 2; task = 'Complex Task'; errass = false; neq = 2; lenwrk = 32*neq; y1keep(1,:) = [ystart]; y2keep(1,:) = zeros(1,2); y3keep(1,:) = zeros(1,2); tkeep(1,:) = zeros(1,2); fprintf('d02pd example program results\n\n'); % We run through the calculation twice with two tolerance values; we % output and plot the results for both of them. tol = [1e-04; 1e-05]; for itol = 1:2 % d02pv is a setup routine to be called prior to d02pd. [work, ifail] = d02pv(tstart, ystart, tend, tol(itol), thres, ... int64(method), task, errass, int64(lenwrk), 'neq', ... int64(neq)); if ifail ~= 0 % Illegal input. Print message and exit. error('Warning: d02pv returned with ifail = %1d ',ifail); end fprintf('Calculation with tol = %1.3e\n\n',tol(itol)); disp(' t y1 y1'''); fprintf('%1.3f %8.4f %8.4f\n', tstart, ystart(1), ystart(2)); % Initialize the storage counter, and the current t value before % looping over t. ncall = 0; tnow = tstart; while (tnow < tend) [tnow, ynow, ypnow, work, ifail] = d02pd(@f, ... int64(neq), work); if ifail ~= 0 % Illegal input. Print message and exit. error('Warning: d02pd returned with ifail = %1d ',... ifail); end % Output results, then store them for plotting. fprintf('%1.3f %8.4f %8.4f\n', tnow, ynow(1), ynow(2)); ncall = ncall + 1; tkeep(ncall,itol) = tnow; if (itol == 2) y1keep(ncall,:) = [ynow(1) ynow(2)]; y2keep(ncall,:) = [ynow(1)-sin(tnow) ynow(2)-cos(tnow)]; else y3keep(ncall,:) = [ynow(1)-sin(tnow) ynow(2)-cos(tnow)]; end end % d02py is a diagnostic routine. [totfcn, stpcst, waste, stpsok, hnext, ifail] = d02py; fprintf(['\nCost of the integration in evaluations of F is ' ... '%1.0f\n\n'], totfcn); end % Plot results. fig = figure('Number', 'off'); display_plot(tkeep, y1keep, y2keep, y3keep) 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,y3keep) % 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. % Calculate spline through data values. xx = tkeep(1,2):0.01:tkeep(length(tkeep(:,2)),2); yy1 = spline(tkeep(:,2),y1keep(:,1),xx); yy2 = spline(tkeep(:,2),y1keep(:,2),xx); % Plot two of the curves and then add the other three. [haxes, hline1, hline2] = plotyy(xx, yy2, tkeep(1:length(y3keep),1), ... abs(y3keep(:,1)), 'plot', 'semilogy'); hold on plot(xx,yy1,':'); hold on plot(tkeep(:,2),y1keep(:,1),'+'); hold on plot(tkeep(:,2),y1keep(:,2),'x'); % Select the required y axis, and plot the sixth curve. axes(haxes(2)); hold on; hline3 = plot(tkeep(:,2), abs(y2keep(:,1))); % 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 1.2]); set(haxes(2), 'YLim', [1e-9 1e-4]); set(haxes(2), 'YMinorTick', 'on'); set(haxes(2), 'YTick', [1e-9 1e-8 1e-7 1e-6 1e-5 1e-4]); 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 Step-by-step RK ',... 'Medium-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.0001)','y-error (tol= 0.00001)','cos(x)',... 'sin(x)','y-solution','y''-solution','Location','Best'); % Set some features of the lines. set(hline1, 'Linewidth', 0.5,'LineStyle','--'); set(hline2, 'Linewidth', 0.5, 'Marker', 's','LineStyle','-.'); set(hline3, 'Linewidth', 0.5, 'Marker', '*','LineStyle',':');
d02pd example program results Calculation with tol = 1.000e-04 t y1 y1' 0.000 0.0000 1.0000 0.785 0.7071 0.7071 1.519 0.9987 0.0513 2.282 0.7573 -0.6531 2.911 0.2285 -0.9735 3.706 -0.5348 -0.8450 4.364 -0.9399 -0.3414 5.320 -0.8209 0.5710 5.802 -0.4631 0.8863 6.283 0.0000 1.0000 Cost of the integration in evaluations of F is 78 Calculation with tol = 1.000e-05 t y1 y1' 0.000 0.0000 1.0000 0.393 0.3827 0.9239 0.785 0.7071 0.7071 1.416 0.9881 0.1538 1.870 0.9557 -0.2943 2.204 0.8062 -0.5916 2.761 0.3711 -0.9286 3.230 -0.0880 -0.9961 3.587 -0.4304 -0.9026 4.022 -0.7710 -0.6368 4.641 -0.9974 -0.0717 5.152 -0.9049 0.4256 5.521 -0.6903 0.7235 5.902 -0.3718 0.9283 6.283 0.0000 1.0000 Cost of the integration in evaluations of F is 118