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′ = 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_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 method = 3${\mathbf{method}}=3$ or -3$-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:     n – int64int32nag_int scalar
n$n$, the number of ordinary differential equations in the system to be solved by the integration function.
Constraint: n1${\mathbf{n}}\ge 1$.
2:     twant – double scalar
t$t$, the value of the independent variable where a solution is desired.
3:     ideriv – int64int32nag_int scalar
Determines whether the solution and/or its first derivative are to be computed
ideriv = 0${\mathbf{ideriv}}=0$
compute approximate solution.
ideriv = 1${\mathbf{ideriv}}=1$
compute approximate first derivative.
ideriv = 2${\mathbf{ideriv}}=2$
compute approximate solution and first derivative.
Constraint: ideriv = 0${\mathbf{ideriv}}=0$, 1$1$ or 2$2$.
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 {\mathbf{n}}$.
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_rkts_onestep (d02pf).
[yp, user] = f(t, n, y, user)

Input Parameters

1:     t – double scalar
t$t$, the current value of the independent variable.
2:     n – int64int32nag_int scalar
n$\mathit{n}$, the number of ordinary differential equations in the system to be solved.
3:     y(n) – 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}$.
4:     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:     yp(n) – double array
The values of fi${f}_{\mathit{i}}$, for i = 1,2,,n$\mathit{i}=1,2,\dots ,\mathit{n}$.
2:     user – Any MATLAB object
6:     wcomm(lwcomm) – double array
This array stores information that can be utilized on subsequent calls to nag_ode_ivp_rkts_interp (d02ps).
7:     iwsav(130$130$) – int64int32nag_int array
8:     rwsav(32 × n + 350$32×{\mathbf{n}}+350$) – 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:     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):
• method = 1${\mathbf{method}}=1$ or -1$-1$ then lwcomm must be at least 1$1$.
• method = 2${\mathbf{method}}=2$ or -2$-2$ then lwcomm must be at least n + max (n,5 × nwant)${\mathbf{n}}+\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left({\mathbf{n}},5×{\mathbf{nwant}}\right)$.
• method = 3${\mathbf{method}}=3$ or -3$-3$ then wcomm and lwcomm are not referenced.
2:     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.

iuser ruser

Output Parameters

1:     ywant(nwant) – double array
An approximation to the first nwant components of the solution at twant if ideriv = 0${\mathbf{ideriv}}=0$ or 2$2$. Otherwise ywant is not defined.
2:     ypwant(nwant) – double array
An approximation to the first nwant components of the first derivative at twant if ideriv = 1${\mathbf{ideriv}}=1$ or 2$2$. Otherwise ypwant is not defined.
3:     wcomm(lwcomm) – double array
Communication array, used to store information between calls to nag_ode_ivp_rkts_interp (d02ps).
4:     user – Any MATLAB object
5:     iwsav(130$130$) – int64int32nag_int array
6:     rwsav(32 × n + 350$32×{\mathbf{n}}+350$) – 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:     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$
Constraint: 1nwantn$1\le {\mathbf{nwant}}\le {\mathbf{n}}$.
Constraint: for method = -1${\mathbf{method}}=-1$ or 1$1$, lwcomm1${\mathbf{lwcomm}}\ge 1$.
Constraint: for method = -2${\mathbf{method}}=-2$ or 2$2$, lwcommn + max (n,5 × nwant)${\mathbf{lwcomm}}\ge {\mathbf{n}}+\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left({\mathbf{n}},5×{\mathbf{nwant}}\right)$.
Constraint: ideriv = 0${\mathbf{ideriv}}=0$, 1$1$ or 2$2$.
method = -3${\mathbf{method}}=-3$ or 3$3$ in setup, but interpolation is not available for this method. Either use method = -2${\mathbf{method}}=-2$ or 2$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, n = _${\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.

Accuracy

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

None.

Example

```function nag_ode_ivp_rkts_interp_example
% 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 = 16;
wcomm = zeros(n+5*nwant, 1);
twant = zeros(npts+1, 1);
ywant = zeros(npts, nwant);
ypwant = zeros(npts, nwant);

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

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

tol = tol*0.1;

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

fprintf('\nCalculation with TOL = %8.1e\n\n', tol);
fprintf('    t         y1         y1''\n');
fprintf(' %6.3f   %8.4f   %8.4f\n', tstart, yinit);

% 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);
tnow = tstart;
j=1;

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

% Integrate one step to tnow.
[tnow, ynow, ypnow, user, iwsav, rwsav, ifail] = ...
nag_ode_ivp_rkts_onestep(@f, n, iwsav, rwsav);

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

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

% Get integration statistics.
[fevals, stepcost, waste, stepsok, hnext, iwsav, ifail] = ...
nag_ode_ivp_rkts_diag(iwsav, rwsav);
fprintf('Cost of the integration in evaluations of f is %d\n', fevals);

end

% Plot results
fig = figure('Number', 'off');
title('Simple Sine Solution, tol=1e-4');
hold on;
%axis([0 10 -1.2 1.2]);
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');

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

Calculation with TOL =  1.0e-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 152

Calculation with TOL =  1.0e-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 231

```
```function d02ps_example
% 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 = 16;
wcomm = zeros(n+5*nwant, 1);
twant = zeros(npts+1, 1);
ywant = zeros(npts, nwant);
ypwant = zeros(npts, nwant);

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

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

fprintf('\nCalculation with TOL = %8.1e\n\n', tol);
fprintf('    t         y1         y1''\n');
fprintf(' %6.3f   %8.4f   %8.4f\n', tstart, yinit);

% 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);
tnow = tstart;
j=1;

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

% Integrate one step to tnow.
[tnow, ynow, ypnow, user, iwsav, rwsav, ifail] = ...
d02pf(@f, n, iwsav, rwsav);

% Interpolate at required additional points up to tnow
while twant(j+1) < tnow
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);
fprintf(' %6.3f   %8.4f   %8.4f\n', twant(j), ywant(j, 1), ypwant(j, 1));

% 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('Cost of the integration in evaluations of f is %d\n', fevals);

end

% Plot results
fig = figure('Number', 'off');
title('Simple Sine Solution, tol=1e-4');
hold on;
%axis([0 10 -1.2 1.2]);
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');

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

Calculation with TOL =  1.0e-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 152

Calculation with TOL =  1.0e-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 231

```