Integer type:  int32  int64  nag_int  show int32  show int32  show int64  show int64  show nag_int  show nag_int

Chapter Contents
Chapter Introduction
NAG Toolbox

# NAG Toolbox: nag_ode_ivp_rkts_interp (d02ps)

## Purpose

nag_ode_ivp_rkts_interp (d02ps) computes the solution of a system of ordinary differential equations using interpolation anywhere on an integration step taken by nag_ode_ivp_rkts_onestep (d02pf).

## Syntax

[ywant, ypwant, wcomm, user, iwsav, rwsav, ifail] = d02ps(n, twant, ideriv, nwant, f, wcomm, iwsav, rwsav, 'lwcomm', lwcomm, 'user', user)
[ywant, ypwant, wcomm, user, iwsav, rwsav, ifail] = nag_ode_ivp_rkts_interp(n, twant, ideriv, nwant, f, wcomm, iwsav, rwsav, 'lwcomm', lwcomm, 'user', user)

## Description

nag_ode_ivp_rkts_interp (d02ps) and its associated functions (nag_ode_ivp_rkts_onestep (d02pf), nag_ode_ivp_rkts_setup (d02pq), nag_ode_ivp_rkts_reset_tend (d02pr), nag_ode_ivp_rkts_diag (d02pt) and nag_ode_ivp_rkts_errass (d02pu)) solve the initial value problem for a first-order system of ordinary differential equations. The functions, based on Runge–Kutta methods and derived from RKSUITE (see Brankin et al. (1991)), integrate
 $y′=ft,y given yt0=y0$
where $y$ is the vector of $\mathit{n}$ solution components and $t$ is the independent variable.
nag_ode_ivp_rkts_onestep (d02pf) computes the solution at the end of an integration step. Using the information computed on that step nag_ode_ivp_rkts_interp (d02ps) computes the solution by interpolation at any point on that step. It cannot be used if ${\mathbf{method}}=3$ or $-3$ was specified in the call to setup function nag_ode_ivp_rkts_setup (d02pq).

## References

Brankin R W, Gladwell I and Shampine L F (1991) RKSUITE: A suite of Runge–Kutta codes for the initial value problems for ODEs SoftReport 91-S1 Southern Methodist University

## Parameters

### Compulsory Input Parameters

1:     $\mathrm{n}$int64int32nag_int scalar
$n$, the number of ordinary differential equations in the system to be solved by the integration function.
Constraint: ${\mathbf{n}}\ge 1$.
2:     $\mathrm{twant}$ – double scalar
$t$, the value of the independent variable where a solution is desired.
3:     $\mathrm{ideriv}$int64int32nag_int scalar
Determines whether the solution and/or its first derivative are to be computed
${\mathbf{ideriv}}=0$
compute approximate solution.
${\mathbf{ideriv}}=1$
compute approximate first derivative.
${\mathbf{ideriv}}=2$
compute approximate solution and first derivative.
Constraint: ${\mathbf{ideriv}}=0$, $1$ or $2$.
4:     $\mathrm{nwant}$int64int32nag_int scalar
The number of components of the solution to be computed. The first nwant components are evaluated.
Constraint: $1\le {\mathbf{nwant}}\le {\mathbf{n}}$.
5:     $\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}$. It must be the same procedure as supplied to nag_ode_ivp_rkts_onestep (d02pf).
[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_interp (d02ps) with the object supplied to nag_ode_ivp_rkts_interp (d02ps).

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
6:     $\mathrm{wcomm}\left({\mathbf{lwcomm}}\right)$ – double array
This array stores information that can be utilized on subsequent calls to nag_ode_ivp_rkts_interp (d02ps).
7:     $\mathrm{iwsav}\left(130\right)$int64int32nag_int array
8:     $\mathrm{rwsav}\left(32×{\mathbf{n}}+350\right)$ – double array
These must be the same arrays supplied in a previous call nag_ode_ivp_rkts_onestep (d02pf). They must remain unchanged between calls.

### Optional Input Parameters

1:     $\mathrm{lwcomm}$int64int32nag_int scalar
Default: the dimension of the array wcomm.
Length of wcomm.
If in a previous call to nag_ode_ivp_rkts_setup (d02pq):
• ${\mathbf{method}}=1$ or $-1$ then lwcomm must be at least $1$.
• ${\mathbf{method}}=2$ or $-2$ then lwcomm must be at least ${\mathbf{n}}+\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left({\mathbf{n}},5×{\mathbf{nwant}}\right)$.
• ${\mathbf{method}}=3$ or $-3$ then wcomm and lwcomm are not referenced.
2:     $\mathrm{user}$ – Any MATLAB object
user is not used by nag_ode_ivp_rkts_interp (d02ps), 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{ywant}\left({\mathbf{nwant}}\right)$ – double array
An approximation to the first nwant components of the solution at twant if ${\mathbf{ideriv}}=0$ or $2$. Otherwise ywant is not defined.
2:     $\mathrm{ypwant}\left({\mathbf{nwant}}\right)$ – double array
An approximation to the first nwant components of the first derivative at twant if ${\mathbf{ideriv}}=1$ or $2$. Otherwise ypwant is not defined.
3:     $\mathrm{wcomm}\left({\mathbf{lwcomm}}\right)$ – double array
Communication array, used to store information between calls to nag_ode_ivp_rkts_interp (d02ps).
4:     $\mathrm{user}$ – Any MATLAB object
5:     $\mathrm{iwsav}\left(130\right)$int64int32nag_int array
6:     $\mathrm{rwsav}\left(32×{\mathbf{n}}+350\right)$ – double array
Information about the integration for use on subsequent calls to nag_ode_ivp_rkts_onestep (d02pf), nag_ode_ivp_rkts_interp (d02ps) 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

Errors or warnings detected by the function:
${\mathbf{ifail}}=1$
Constraint: $1\le {\mathbf{nwant}}\le {\mathbf{n}}$.
Constraint: for ${\mathbf{method}}=-1$ or $1$, ${\mathbf{lwcomm}}\ge 1$.
Constraint: for ${\mathbf{method}}=-2$ or $2$, ${\mathbf{lwcomm}}\ge {\mathbf{n}}+\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left({\mathbf{n}},5×{\mathbf{nwant}}\right)$.
Constraint: ${\mathbf{ideriv}}=0$, $1$ or $2$.
${\mathbf{method}}=-3$ or $3$ in setup, but interpolation is not available for this method. Either use ${\mathbf{method}}=-2$ or $2$ in setup or use reset function to force the integrator to step to particular points.
On entry, a previous call to the setup function has not been made or the communication arrays have become corrupted, or a catastrophic error has already been detected elsewhere.
You cannot continue integrating the problem.
On entry, ${\mathbf{n}}=_$, but the value passed to the setup function was .
You cannot call this function after the integrator has returned an error.
You cannot call this function before you have called the step integrator.
You cannot call this function when you have specified, in the setup function, that the range integrator will be used.
${\mathbf{ifail}}=-99$
${\mathbf{ifail}}=-399$
Your licence key may have expired or may not have been installed correctly.
${\mathbf{ifail}}=-999$
Dynamic memory allocation failed.

## Accuracy

The computed values will be of a similar accuracy to that computed by nag_ode_ivp_rkts_onestep (d02pf).

None.

## Example

This example solves the equation
 $y′′ = -y , y0=0, y′0=1$
reposed as
 $y1′ = y2$
 $y2′ = -y1$
over the range $\left[0,2\pi \right]$ with initial conditions ${y}_{1}=0.0$ and ${y}_{2}=1.0$. Relative error control is used with threshold values of $\text{1.0e−8}$ for each solution component. nag_ode_ivp_rkts_onestep (d02pf) is used to integrate the problem one step at a time and nag_ode_ivp_rkts_interp (d02ps) is used to compute the first component of the solution and its derivative at intervals of length $\pi /8$ across the range whenever these points lie in one of those integration steps. A low order Runge–Kutta method (${\mathbf{method}}=-1$) is also used with tolerances ${\mathbf{tol}}=\text{1.0e−4}$ and ${\mathbf{tol}}=\text{1.0e−5}$ in turn so that solutions may be compared.
```function d02ps_example

fprintf('d02ps 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);
nwant  = int64(1);
tol    =  1.0E-2;
npts   = 32;
wcomm  = zeros(n+5*nwant, 1);
twant  = zeros(npts+1, 1);
ywant  = zeros(npts, nwant);
ypwant = zeros(npts, nwant);

% Set control for printing solution
tinc   = (tend-tstart)/npts;

% 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);

% Set up first point at which solution is required
twant(1) = tstart;
twant(2) = tstart + tinc;
ywant(1,1) = yinit(1);
ypwant(1,1) = yinit(2);
t = tstart;
j=1;

% Integrate by steps until tend is reached or error is encountered.
while t < tend

% Integrate one step from t, updating t.
[t, y, yp, user, iwsav, rwsav, ifail] = d02pf(@f, n, iwsav, rwsav);

% Interpolate at required additional points up to t
while twant(j+1) < t
j=j+1;
% Interpolate and print solution at t = twant.
ideriv = int64(2);
[ywant(j, :), ypwant(j, :), wcomm, user, iwsav, rwsav, ifail] = ...
d02ps(n, twant(j), ideriv, nwant, @f, wcomm, iwsav, rwsav);

% Set next required solution point
twant(j+1) = twant(j) + tinc;
end
end

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

end

% Plot results
fig1 = figure;
title('Simple Sine Solution, tol=1e-4');
hold on;
xlabel('t');
ylabel('Solution');
plot(twant(1:npts), ywant(:, 1), '-xr');
text(3, 0.25, 'Solution', 'Color', 'r');
plot(twant(1:npts), ypwant(:, 1), '-xg');
text(1.45, 0.25, 'Derivative', 'Color', 'g');
hold off

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

Calculation with TOL =  1.0e-03

Number of evaluations of f = 152

Calculation with TOL =  1.0e-04

Number of evaluations of f = 231
```