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_withdraw_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_withdraw_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′=ft,y given yt0=y0$
where $y$ is the vector of $\mathit{n}$ solution components and $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 ${\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:     $\mathrm{neq}$int64int32nag_int scalar
$n$, the number of ordinary differential equations in the system to be solved by the integration function.
Constraint: ${\mathbf{neq}}\ge 1$.
2:     $\mathrm{twant}$ – double scalar
$t$, the value of the independent variable where a solution is desired.
3:     $\mathrm{reqest}$ – string (length ≥ 1)
Determines whether the solution and/or its first derivative are to be computed.
${\mathbf{reqest}}=\text{'S'}$
Compute the approximate solution only.
${\mathbf{reqest}}=\text{'D'}$
Compute the approximate first derivative of the solution only.
${\mathbf{reqest}}=\text{'B'}$
Compute both the approximate solution and its first derivative.
Constraint: ${\mathbf{reqest}}=\text{'S'}$, $\text{'D'}$ or $\text{'B'}$.
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 \mathit{n}$, where $\mathit{n}$ is specified by neq in the prior call to nag_ode_ivp_rk_setup (d02pv).
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_rk_onestep (d02pd).
[yp] = f(t, y)

Input Parameters

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

Output Parameters

1:     $\mathrm{yp}\left(:\right)$ – double array
The values of ${f}_{\mathit{i}}$, for $\mathit{i}=1,2,\dots ,\mathit{n}$.
6:     $\mathrm{work}\left(:\right)$ – double array
The dimension of the array work must be at least ${\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:     $\mathrm{wrkint}\left({\mathbf{lenint}}\right)$ – double array
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:     $\mathrm{lenint}$int64int32nag_int scalar
Default: the dimension of the array wrkint.
The dimension of the array wrkint.
Constraints:
• ${\mathbf{lenint}}\ge 1$ if ${\mathbf{method}}=1$ in the prior call to nag_ode_ivp_rk_setup (d02pv);
• ${\mathbf{lenint}}\ge \mathit{n}+5×{\mathbf{nwant}}$ if ${\mathbf{method}}=2$ and $\mathit{n}$ is specified by neq in the prior call to nag_ode_ivp_rk_setup (d02pv).

### Output Parameters

1:     $\mathrm{ywant}\left(:\right)$ – double array
The dimension of the array ywant will be ${\mathbf{nwant}}$ if ${\mathbf{reqest}}=\text{'S'}$ or $\text{'B'}$ and $1$ otherwise
An approximation to the first nwant components of the solution at twant if ${\mathbf{reqest}}=\text{'S'}$ or $\text{'B'}$. Otherwise ywant is not defined.
2:     $\mathrm{ypwant}\left(:\right)$ – double array
The dimension of the array ypwant will be ${\mathbf{nwant}}$ if ${\mathbf{reqest}}=\text{'D'}$ or $\text{'B'}$ and $1$ otherwise
An approximation to the first nwant components of the first derivative at twant if ${\mathbf{reqest}}=\text{'D'}$ or $\text{'B'}$. Otherwise ypwant is not defined.
3:     $\mathrm{work}\left(:\right)$ – double array
The dimension of the array work will be ${\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:     $\mathrm{wrkint}\left({\mathbf{lenint}}\right)$ – double array
The contents are modified.
5:     $\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$
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 ${\mathbf{method}}=3$. You cannot continue integrating the problem.
${\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_rk_onestep (d02pd).

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_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 $\pi /8$ across the range whenever these points lie in one of those integration steps. A moderate order Runge–Kutta method (${\mathbf{method}}=2$) is also used with tolerances ${\mathbf{tol}}=\text{1.0e−3}$ and ${\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 argument nwant.
```function d02px_example

fprintf('d02px 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);
errass = false;
lenwrk = int64(32*n);
reqest = 'Both';
nwant  = int64(1);
tol    =  1.0E-2;
npts   = 32;
work   = zeros(n+5*nwant, 1);
twant  = zeros(npts+1, 1);
ywant  = zeros(npts, nwant);
ypwant = zeros(npts, nwant);
lenint = n+5*nwant;
wrkint = zeros(lenint,1);

% 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
[work, ifail] = d02pv(tstart, yinit, tend, tol, thresh, method, task, ...
errass, lenwrk);

% 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, work, ifail] = d02pd(@f, n, work);

% Interpolate at required additional points up to t
while twant(j+1) < t
j=j+1;
% Interpolate and print solution at t = twant.
[ywant(j,:), ypwant(j,:), work, wrkint, ifail] = d02px(n, twant(j), ...
reqest, nwant, @f, work, wrkint);

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

% Get integration statistics.
[fevals, stepcost, waste, stepsok, hnext, fail] = d02py;
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] = f(t, y)
% Evaluate derivative vector.
yp = zeros(2, 1);
yp(1) =  y(2);
yp(2) = -y(1);
```
```d02px example results

Calculation with TOL =  1.0e-03

Number of evaluations of f = 68

Calculation with TOL =  1.0e-04

Number of evaluations of f = 105
```