NAG Toolbox: nag_ode_ivp_rkts_onestep (d02pf)

Purpose

Syntax

Description

References

Parameters

Compulsory Input Parameters

1:     $\mathrm{f}$ – function handle or string containing name of m-file

f must evaluate the functions $f_i$ (that is the first derivatives $y_i^{\prime }$) for given values of the arguments $t$, $y_i$.

[yp, user] = f(t, n, y, user)

Input Parameters

1:     $\mathrm{t}$ – double scalar

$t$, the current value of the independent variable.

2:     $\mathrm{n}$ – int64int32nag_int scalar

$\mathit{n}$, the number of ordinary differential equations in the system to be solved.

3:     $\mathrm{y}\left(\mathbf{n}\right)$ – double array

The current values of the dependent variables, $y_{\mathit{i}}$, for $\mathit{i}=1,2,\dots ,\mathit{n}$.

4:     $\mathrm{user}$ – Any MATLAB object

f is called from nag_ode_ivp_rkts_onestep (d02pf) with the object supplied to nag_ode_ivp_rkts_onestep (d02pf).

Output Parameters

1:     $\mathrm{yp}\left(\mathbf{n}\right)$ – double array

The values of $f_{\mathit{i}}$, for $\mathit{i}=1,2,\dots ,\mathit{n}$.

2:     $\mathrm{user}$ – Any MATLAB object

2:     $\mathrm{n}$ – int64int32nag_int scalar

$n$, the number of ordinary differential equations in the system to be solved.

Constraint: $\mathbf{n}\ge 1$.

3:     $\mathrm{iwsav}\left(130\right)$ – int64int32nag_int array

4:     $\mathrm{rwsav}\left(32\times \mathbf{n}+350\right)$ – double array

These must be the same arrays supplied in a previous call to nag_ode_ivp_rkts_setup (d02pq). They must remain unchanged between calls.

Optional Input Parameters

1:     $\mathrm{user}$ – Any MATLAB object

user is not used by nag_ode_ivp_rkts_onestep (d02pf), but is passed to f. Note that for large objects it may be more efficient to use a global variable which is accessible from the m-files than to use user.

Output Parameters

1:     $\mathrm{tnow}$ – double scalar

$t$, the value of the independent variable at which a solution has been computed.

2:     $\mathrm{ynow}\left(\mathbf{n}\right)$ – double array

An approximation to the solution at tnow. The local error of the step to tnow was no greater than permitted by the specified tolerances (see nag_ode_ivp_rkts_setup (d02pq)).

3:     $\mathrm{ypnow}\left(\mathbf{n}\right)$ – double array

An approximation to the first derivative of the solution at tnow.

4:     $\mathrm{user}$ – Any MATLAB object

5:     $\mathrm{iwsav}\left(130\right)$ – int64int32nag_int array

6:     $\mathrm{rwsav}\left(32\times \mathbf{n}+350\right)$ – double array

Information about the integration for use on subsequent calls to nag_ode_ivp_rkts_onestep (d02pf) or other associated functions.

7:     $\mathrm{ifail}$ – int64int32nag_int scalar

$\mathbf{ifail}=\mathbf{0}$ unless the function detects an error (see Error Indicators and Warnings).

Error Indicators and Warnings

Cases prefixed with W are classified as warnings and do not generate an error of type NAG:error_n. See nag_issue_warnings.

   $\mathbf{ifail}=1$

A call to this function cannot be made after it has returned an error.
The setup function must be called to start another problem.

On entry, $\mathbf{n}=\_$, but the value passed to the setup function was .

On entry, the communication arrays have become corrupted, or a catastrophic error has already been detected elsewhere. You cannot continue integrating the problem.

tend, as specified in the setup function, has already been reached.
To start a new problem, you will need to call the setup function.
To continue integration beyond tend then
nag_ode_ivp_rkts_reset_tend (d02pr) must first be called to reset tend to a new end value.

W  $\mathbf{ifail}=2$

More than $100$ output points have been obtained by integrating to tend (as specified in the setup function). They have been so clustered that it would probably be (much) more efficient to use the interpolation function (if $\left|\mathbf{method}\right|=3$, switch to $\left|\mathbf{method}\right|=2$ at setup).
However, you can continue integrating the problem.

W  $\mathbf{ifail}=3$

Approximately $\_$ function evaluations have been used to compute the solution since the integration started or since this message was last printed.
However, you can continue integrating the problem.

W  $\mathbf{ifail}=4$

Approximately $\_$ function evaluations have been used to compute the solution since the integration started or since this message was last printed. Your problem has been diagnosed as stiff. If the situation persists, it will cost roughly $\_$ times as much to reach tend (setup) as it has cost to reach the current time. You should probably call functions intended for stiff problems. However, you can continue integrating the problem.

W  $\mathbf{ifail}=5$

In order to satisfy your error requirements the solver has to use a step size of $\_$ at the current time, $\_$. This step size is too small for the machine precision, and is smaller than $\_$.

W  $\mathbf{ifail}=6$

The global error assessment algorithm failed at start of integration.
The integration is being terminated.

The global error assessment may not be reliable for times beyond $\_$.
The integration is being terminated.

   $\mathbf{ifail}=-99$

An unexpected error has been triggered by this routine. Please contact NAG.

   $\mathbf{ifail}=-399$

Your licence key may have expired or may not have been installed correctly.

   $\mathbf{ifail}=-999$

Dynamic memory allocation failed.

Accuracy

Further Comments

Example

Open in the MATLAB editor: d02pf_example

function d02pf_example


fprintf('d02pf example results\n\n');

% Set initial conditions and input
method = int64(2);
tstart = 0;
tend = 2*pi;
yinit = [0;1];
hstart = 0;
thresh = [1e-08; 1e-08];
n = int64(2);
tol0 =  1.0E-4;
ynow = zeros(20, n);
tnow = zeros(20, 1);
err1 = zeros(20, 2);
err2 = zeros(20, 2);

% Set control for output
tol = 10.0*tol0;

% Run through the calculation twice with two tolerance values
for i = 1:2

  tol = tol*0.1;

  % Call setup function
  [iwsav, rwsav, ifail] = d02pq(tstart, tend, yinit, tol, thresh, method);


  tnow(1) = tstart;
  ynow(1,:) = yinit;
  j=1;
  while tnow(j) < tend
    j=j+1;
    [tnow(j), ynow(j, :), ypnow, user, iwsav, rwsav, ifail] = ...
      d02pf(@f, n, iwsav, rwsav);

    err1(j, i) =  ynow(j, 1)-sin(tnow(j));
    err2(j, i) =  ynow(j, 2)-cos(tnow(j));

  end

  fprintf('\nCalculation with TOL = %8.1e:\n\n', tol);
  [fevals, stepcost, waste, stepsok, hnext, iwsav, ifail] = d02pt(iwsav, rwsav);
  fprintf('  Number of evaluations of f = %d\n', fevals);

  % Store the t values for use in plotting errors
  if i == 1
    tnow1 = tnow;
  end

  % Store value of j
  npts(i) = j;
end

% Plot results
fig1 = figure;
title({['First-order ODEs solution by single stepping'],...
          ['Medium-order Runge-Kutta Method, Two Tolerances']});
hold on;
axis([0 10 -1.2 1.2]);
xlabel('t');
ylabel('Solution (y, y'')');
plot(tnow(1:npts(2)), ynow(1:npts(2), 1), '-xr');
text(ceil(tnow(npts(2))), ynow(npts(2), 1), 'y', 'Color', 'r');
plot(tnow(1:npts(2)), ynow(1:npts(2), 2), '-xg');
text(ceil(tnow(npts(2))), ynow(npts(2), 2), 'y''', 'Color', 'g');
% Plot errors with a different (log) scale
ax1 = gca;
ax2 = axes('Position',get(ax1,'Position'),...
           'XAxisLocation','bottom','YAxisLocation','right',...
           'YScale','log','Color','none','XColor','k','YColor','k');
hold on;
axis([0 10 1e-9 1e-4]);
ylabel('abs(Error)');
plot(ax2, tnow1(1:npts(1)), abs(err1(1:npts(1), 1)), '-*b');
text(ceil(tnow1(npts(1))), err1(npts(1), 1) - 1e-5, ...
     'y-error (tol=0.001)', 'Color', 'b');
plot(ax2, tnow, abs(err1(:, 2)), '-sm');
text(ceil(tnow(npts(2))), err1(npts(2), 2), 'y-error (tol=0.0001)', ...
     'Color', 'm');
hold off;



function [yp, user] = f(t, n, y, user)
  yp = [y(2); -y(1)];

d02pf example results


Calculation with TOL =  1.0e-04:

  Number of evaluations of f = 204

Calculation with TOL =  1.0e-05:

  Number of evaluations of f = 314