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_rk_interp (d02px)

## Purpose

nag_ode_ivp_rk_interp (d02px) computes the solution of a system of ordinary differential equations using interpolation anywhere on an integration step taken by nag_ode_ivp_rk_onestep (d02pd).
Note: this function is scheduled to be withdrawn, please see d02px in Advice on Replacement Calls for Withdrawn/Superseded Routines..

## Syntax

[ywant, ypwant, work, wrkint, ifail] = d02px(neq, twant, reqest, nwant, f, work, wrkint, 'lenint', lenint)
[ywant, ypwant, work, wrkint, ifail] = nag_ode_ivp_rk_interp(neq, twant, reqest, nwant, f, work, wrkint, 'lenint', lenint)

## Description

nag_ode_ivp_rk_interp (d02px) and its associated functions (nag_ode_ivp_rk_onestep (d02pd), nag_ode_ivp_rk_setup (d02pv), nag_ode_ivp_rk_reset_tend (d02pw), nag_ode_ivp_rk_diag (d02py) and nag_ode_ivp_rk_errass (d02pz)) 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′ = f(t,y)  given  y(t0) = y0 $y′=f(t,y) given y(t0)=y0$
where y$y$ is the vector of n$\mathit{n}$ solution components and t$t$ is the independent variable.
nag_ode_ivp_rk_onestep (d02pd) computes the solution at the end of an integration step. Using the information computed on that step nag_ode_ivp_rk_interp (d02px) computes the solution by interpolation at any point on that step. It cannot be used if method = 3${\mathbf{method}}=3$ was specified in the call to setup function nag_ode_ivp_rk_setup (d02pv).

## 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:     neq – int64int32nag_int scalar
n$n$, the number of ordinary differential equations in the system to be solved by the integration function.
Constraint: neq1${\mathbf{neq}}\ge 1$.
2:     twant – double scalar
t$t$, the value of the independent variable where a solution is desired.
3:     reqest – string (length ≥ 1)
Determines whether the solution and/or its first derivative are to be computed.
reqest = 'S'${\mathbf{reqest}}=\text{'S'}$
Compute the approximate solution only.
reqest = 'D'${\mathbf{reqest}}=\text{'D'}$
Compute the approximate first derivative of the solution only.
reqest = 'B'${\mathbf{reqest}}=\text{'B'}$
Compute both the approximate solution and its first derivative.
Constraint: reqest = 'S'${\mathbf{reqest}}=\text{'S'}$, 'D'$\text{'D'}$ or 'B'$\text{'B'}$.
4:     nwant – int64int32nag_int scalar
The number of components of the solution to be computed. The first nwant components are evaluated.
Constraint: 1nwantn$1\le {\mathbf{nwant}}\le \mathit{n}$, where n$\mathit{n}$ is specified by neq in the prior call to nag_ode_ivp_rk_setup (d02pv).
5:     f – function handle or string containing name of m-file
f must evaluate the functions fi${f}_{i}$ (that is the first derivatives yi${y}_{i}^{\prime }$) for given values of the arguments t,yi$t,{y}_{i}$. It must be the same procedure as supplied to nag_ode_ivp_rk_onestep (d02pd).
[yp] = f(t, y)

Input Parameters

1:     t – double scalar
t$t$, the current value of the independent variable.
2:     y( : $:$) – double array
The current values of the dependent variables, yi${y}_{\mathit{i}}$, for i = 1,2,,n$\mathit{i}=1,2,\dots ,\mathit{n}$.

Output Parameters

1:     yp( : $:$) – double array
The values of fi${f}_{\mathit{i}}$, for i = 1,2,,n$\mathit{i}=1,2,\dots ,\mathit{n}$.
6:     work( : $:$) – double array
Note: the dimension of the array work must be at least lenwrk${\mathbf{lenwrk}}$ (see nag_ode_ivp_rk_setup (d02pv)).
This must be the same array as supplied to nag_ode_ivp_rk_onestep (d02pd) and must remain unchanged between calls.
7:     wrkint(lenint) – double array
lenint, the dimension of the array, must satisfy the constraint
• lenint1${\mathbf{lenint}}\ge 1$ if method = 1${\mathbf{method}}=1$ in the prior call to nag_ode_ivp_rk_setup (d02pv)
• lenintn + 5 × nwant${\mathbf{lenint}}\ge \mathit{n}+5×{\mathbf{nwant}}$ if method = 2${\mathbf{method}}=2$ and n$\mathit{n}$ is specified by neq in the prior call to nag_ode_ivp_rk_setup (d02pv)
• .
Must be the same array as supplied in previous calls, if any, and must remain unchanged between calls to nag_ode_ivp_rk_interp (d02px).

### Optional Input Parameters

1:     lenint – int64int32nag_int scalar
Default: The dimension of the array wrkint.
The dimension of the array wrkint as declared in the (sub)program from which nag_ode_ivp_rk_interp (d02px) is called.
Constraints:
• lenint1${\mathbf{lenint}}\ge 1$ if method = 1${\mathbf{method}}=1$ in the prior call to nag_ode_ivp_rk_setup (d02pv);
• lenintn + 5 × nwant${\mathbf{lenint}}\ge \mathit{n}+5×{\mathbf{nwant}}$ if method = 2${\mathbf{method}}=2$ and n$\mathit{n}$ is specified by neq in the prior call to nag_ode_ivp_rk_setup (d02pv).

None.

### Output Parameters

1:     ywant( : $:$) – double array
Note: the dimension of the array ywant must be at least nwant${\mathbf{nwant}}$ if reqest = 'S'${\mathbf{reqest}}=\text{'S'}$ or 'B'$\text{'B'}$, and at least 1$1$ otherwise.
An approximation to the first nwant components of the solution at twant if reqest = 'S'${\mathbf{reqest}}=\text{'S'}$ or 'B'$\text{'B'}$. Otherwise ywant is not defined.
2:     ypwant( : $:$) – double array
Note: the dimension of the array ypwant must be at least nwant${\mathbf{nwant}}$ if reqest = 'D'${\mathbf{reqest}}=\text{'D'}$ or 'B'$\text{'B'}$, and at least 1$1$ otherwise.
An approximation to the first nwant components of the first derivative at twant if reqest = 'D'${\mathbf{reqest}}=\text{'D'}$ or 'B'$\text{'B'}$. Otherwise ypwant is not defined.
3:     work( : $:$) – double array
Note: the dimension of the array work must be at least lenwrk${\mathbf{lenwrk}}$ (see nag_ode_ivp_rk_setup (d02pv)).
Contains information about the integration for use on subsequent calls to nag_ode_ivp_rk_onestep (d02pd) or other associated functions.
4:     wrkint(lenint) – double array
The contents are modified.
5:     ifail – int64int32nag_int scalar
${\mathrm{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:
ifail = 1${\mathbf{ifail}}=1$
On entry, an invalid input value for nwant or lenint was detected or an invalid call to nag_ode_ivp_rk_interp (d02px) was made, for example without a previous call to the integration function nag_ode_ivp_rk_onestep (d02pd), or after an error return from nag_ode_ivp_rk_onestep (d02pd), or if nag_ode_ivp_rk_onestep (d02pd) was being used with method = 3${\mathbf{method}}=3$. You cannot continue integrating the problem.

## Accuracy

The computed values will be of a similar accuracy to that computed by nag_ode_ivp_rk_onestep (d02pd).

None.

## Example

This example solves the equation
 y ′ ′ = − y ,   y(0) = 0,   y′(0) = 1 $y′′ = -y , y(0)=0, y′(0)=1$
reposed as
 y1 ′ = y2 $y1′ = y2$
 y2 ′ = − y1 $y2′ = -y1$
over the range [0,2π]$\left[0,2\pi \right]$ with initial conditions y1 = 0.0${y}_{1}=0.0$ and y2 = 1.0${y}_{2}=1.0$. Relative error control is used with threshold values of 1.0e−8$\text{1.0e−8}$ for each solution component. nag_ode_ivp_rk_onestep (d02pd) is used to integrate the problem one step at a time and nag_ode_ivp_rk_interp (d02px) is used to compute the first component of the solution and its derivative at intervals of length π / 8$\pi /8$ across the range whenever these points lie in one of those integration steps. A moderate order Runge–Kutta method (method = 2${\mathbf{method}}=2$) is also used with tolerances tol = 1.0e−3${\mathbf{tol}}=\text{1.0e−3}$ and tol = 1.0e−4${\mathbf{tol}}=\text{1.0e−4}$ in turn so that solutions may be compared. The value of π$\pi$ is obtained by using nag_math_pi (x01aa).
Note that the length of work is large enough for any valid combination of input arguments to nag_ode_ivp_rk_setup (d02pv) and the length of wrkint is large enough for any valid value of the parameter nwant.
```function nag_ode_ivp_rk_interp_example
% Initialize variables and arrays.
tstart = 0;
hstart = 0;
ystart = [0; 1];
tend = 2*pi;
thres = [1e-08; 1e-08];
method = 2;
errass = false;
neq = 2;
lenwrk = 32*neq;
twant = pi/8;
reqest = 'Both';
nwant = 1;
wrkint = zeros(7, 1);
npts = 16;
lenint = neq+5*nwant;
tinc = tend/npts;

fprintf('nag_ode_ivp_rk_interp 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-03; 1e-04];
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
% Unsuccessful call.  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, store the current values.
ncall = 1;
y1keep(ncall) = ystart(1);
y2keep(ncall) = ystart(2);
tkeep(ncall) = tstart;

% Initialize the current t values, and start looping over t. tnow
% nmarks the end of the integration step for nag_ode_ivp_rk_onestep; twant is the
% point on that step at which nag_ode_ivp_rk_interp calculates the solution
% by interpolation.
tnow = tstart;
twant = tinc;
while (tnow < tend)

[tnow, ynow, ypnow, work, ifail] = nag_ode_ivp_rk_onestep(@f, ...
int64(neq), work);
if ifail ~= 0
% Unsuccessful call.  Print message and exit.
error('Warning: nag_ode_ivp_rk_onestep returned with ifail = %1d ',...
ifail);
end

while (twant < tnow)

% Call nag_ode_ivp_rk_interp to obtain solution by interpolation on results
% from nag_ode_ivp_rk_onestep.  Note that the function file must be the same
% for both routines.
[ywant, ypwant, work, wrkint, ifail] = nag_ode_ivp_rk_interp(int64(neq), ...
twant, reqest, int64(nwant), @f, work, wrkint,...
'lenint', int64(lenint));

% Output results, then store them for plotting.
fprintf('%1.3f     %8.4f     %8.4f\n', twant, ywant, ypwant);
ncall = ncall+1;
y1keep(ncall) = ywant;
y2keep(ncall) = ypwant;
tkeep(ncall) = twant;

twant = twant + tinc;
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);

% Plot results.
fig = figure('Number', 'off');
display_plot(tkeep, y1keep, y2keep, tol(itol))
end

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, tol)
% 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 the results.
plot(tkeep, y1keep, '-+', tkeep, y2keep, ':x')
title(['Simple Sine Solution, TOL = ',num2str(tol)], titleFmt{:});
% Label the axes.
xlabel('t', labFmt{:});
ylabel('Solution', labFmt{:});
legend('solution','derivative','Location','Best');
```
```
nag_ode_ivp_rk_interp example program results

Calculation with tol = 1.000e-03

t           y1           y1'
0.000       0.0000       1.0000
0.393       0.3827       0.9239
0.785       0.7071       0.7071
1.178       0.9239       0.3826
1.571       1.0000      -0.0001
1.963       0.9238      -0.3828
2.356       0.7070      -0.7073
2.749       0.3825      -0.9240
3.142      -0.0002      -0.9999
3.534      -0.3829      -0.9238
3.927      -0.7072      -0.7069
4.320      -0.9239      -0.3823
4.712      -0.9999       0.0004
5.105      -0.9236       0.3830
5.498      -0.7068       0.7073
5.890      -0.3823       0.9239

Cost of the integration in evaluations of F is 68

Calculation with tol = 1.000e-04

t           y1           y1'
0.000       0.0000       1.0000
0.393       0.3827       0.9239
0.785       0.7071       0.7071
1.178       0.9239       0.3827
1.571       1.0000      -0.0000
1.963       0.9239      -0.3827
2.356       0.7071      -0.7071
2.749       0.3827      -0.9239
3.142      -0.0000      -1.0000
3.534      -0.3827      -0.9239
3.927      -0.7071      -0.7071
4.320      -0.9239      -0.3827
4.712      -1.0000       0.0000
5.105      -0.9238       0.3827
5.498      -0.7071       0.7071
5.890      -0.3826       0.9239

Cost of the integration in evaluations of F is 105

```
```function d02px_example
% Initialize variables and arrays.
tstart = 0;
hstart = 0;
ystart = [0; 1];
tend = 2*pi;
thres = [1e-08; 1e-08];
method = 2;
errass = false;
neq = 2;
lenwrk = 32*neq;
twant = pi/8;
reqest = 'Both';
nwant = 1;
wrkint = zeros(7, 1);
npts = 16;
lenint = neq+5*nwant;
tinc = tend/npts;

fprintf('d02px 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-03; 1e-04];
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
% Unsuccessful call.  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, store the current values.
ncall = 1;
y1keep(ncall) = ystart(1);
y2keep(ncall) = ystart(2);
tkeep(ncall) = tstart;

% Initialize the current t values, and start looping over t. tnow
% nmarks the end of the integration step for d02pd; twant is the
% point on that step at which d02px calculates the solution
% by interpolation.
tnow = tstart;
twant = tinc;
while (tnow < tend)

[tnow, ynow, ypnow, work, ifail] = d02pd(@f, ...
int64(neq), work);
if ifail ~= 0
% Unsuccessful call.  Print message and exit.
error('Warning: d02pd returned with ifail = %1d ',...
ifail);
end

while (twant < tnow)

% Call d02px to obtain solution by interpolation on results
% from d02pd.  Note that the function file must be the same
% for both routines.
[ywant, ypwant, work, wrkint, ifail] = d02px(int64(neq), ...
twant, reqest, int64(nwant), @f, work, wrkint,...
'lenint', int64(lenint));

% Output results, then store them for plotting.
fprintf('%1.3f     %8.4f     %8.4f\n', twant, ywant, ypwant);
ncall = ncall+1;
y1keep(ncall) = ywant;
y2keep(ncall) = ypwant;
tkeep(ncall) = twant;

twant = twant + tinc;
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);

% Plot results.
fig = figure('Number', 'off');
display_plot(tkeep, y1keep, y2keep, tol(itol))
end

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, tol)
% 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 the results.
plot(tkeep, y1keep, '-+', tkeep, y2keep, ':x')
title(['Simple Sine Solution, TOL = ',num2str(tol)], titleFmt{:});
% Label the axes.
xlabel('t', labFmt{:});
ylabel('Solution', labFmt{:});
legend('solution','derivative','Location','Best');
```
```
d02px example program results

Calculation with tol = 1.000e-03

t           y1           y1'
0.000       0.0000       1.0000
0.393       0.3827       0.9239
0.785       0.7071       0.7071
1.178       0.9239       0.3826
1.571       1.0000      -0.0001
1.963       0.9238      -0.3828
2.356       0.7070      -0.7073
2.749       0.3825      -0.9240
3.142      -0.0002      -0.9999
3.534      -0.3829      -0.9238
3.927      -0.7072      -0.7069
4.320      -0.9239      -0.3823
4.712      -0.9999       0.0004
5.105      -0.9236       0.3830
5.498      -0.7068       0.7073
5.890      -0.3823       0.9239

Cost of the integration in evaluations of F is 68

Calculation with tol = 1.000e-04

t           y1           y1'
0.000       0.0000       1.0000
0.393       0.3827       0.9239
0.785       0.7071       0.7071
1.178       0.9239       0.3827
1.571       1.0000      -0.0000
1.963       0.9239      -0.3827
2.356       0.7071      -0.7071
2.749       0.3827      -0.9239
3.142      -0.0000      -1.0000
3.534      -0.3827      -0.9239
3.927      -0.7071      -0.7071
4.320      -0.9239      -0.3827
4.712      -1.0000       0.0000
5.105      -0.9238       0.3827
5.498      -0.7071       0.7071
5.890      -0.3826       0.9239

Cost of the integration in evaluations of F is 105

```