hide long namesshow long names
hide short namesshow short names
Integer type:  int32  int64  nag_int  show int32  show int32  show int64  show int64  show nag_int  show nag_int

PDF version (NAG web site, 64-bit version, 64-bit version)
Chapter Contents
Chapter Introduction
NAG Toolbox

NAG Toolbox: nag_ode_ivp_rk_setup (d02pv)

Purpose

nag_ode_ivp_rk_setup (d02pv) is a setup function which must be called prior to the first call of either of the integration functions nag_ode_ivp_rk_range (d02pc) and nag_ode_ivp_rk_onestep (d02pd).
Note: this function is scheduled to be withdrawn, please see d02pv in Advice on Replacement Calls for Withdrawn/Superseded Routines..

Syntax

[work, ifail] = d02pv(tstart, ystart, tend, tol, thres, method, task, errass, lenwrk, 'neq', neq, 'hstart', hstart)
[work, ifail] = nag_ode_ivp_rk_setup(tstart, ystart, tend, tol, thres, method, task, errass, lenwrk, 'neq', neq, 'hstart', hstart)

Description

nag_ode_ivp_rk_setup (d02pv) and its associated functions (nag_ode_ivp_rk_range (d02pc), nag_ode_ivp_rk_onestep (d02pd), nag_ode_ivp_rk_reset_tend (d02pw), nag_ode_ivp_rk_interp (d02px), 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 yy is the vector of nn solution components and tt is the independent variable.
The integration proceeds by steps from the initial point t0t0 towards the final point tftf. An approximate solution yy is computed at each step. For each component yiyi, for i = 1,2,,ni=1,2,,n, the error made in the step, i.e., the local error, is estimated. The step size is chosen automatically so that the integration will proceed efficiently while keeping this local error estimate smaller than a tolerance that you specify by means of parameters tol and thres.
nag_ode_ivp_rk_range (d02pc) can be used to solve the ‘usual task’, namely integrating the system of differential equations to obtain answers at points you specify. nag_ode_ivp_rk_onestep (d02pd) is used for all more ‘complicated tasks’.
You should consider carefully how you want the local error to be controlled. Essentially the code uses relative local error control, with tol being the desired relative accuracy. For reliable computation, the code must work with approximate solutions that have some correct digits, so there is an upper bound on the value you can specify for tol. It is impossible to compute a numerical solution that is more accurate than the correctly rounded value of the true solution, so you are not allowed to specify tol too small for the precision you are using. The magnitude of the local error in yiyi on any step will not be greater than tol × max (μi,thres(i))tol×max(μi,thresi) where μiμi is an average magnitude of yiyi over the step. If thres(i)thresi is smaller than the current value of μiμi, this is a relative error test and tol indicates how many significant digits you want in yiyi. If thres(i)thresi is larger than the current value of μiμi, this is an absolute error test with tolerance tol × thres(i)tol×thresi. Relative error control is the recommended mode of operation, but pure relative error control, thres(i) = 0.0thresi=0.0, is not permitted. See Section [Further Comments] for further information about error control.
nag_ode_ivp_rk_range (d02pc) and nag_ode_ivp_rk_onestep (d02pd) control local error rather than the true (global) error, the difference between the numerical and true solution. Control of the local error controls the true error indirectly. Roughly speaking, the code produces a solution that satisfies the differential equation with a discrepancy bounded in magnitude by the error tolerance. What this implies about how close the numerical solution is to the true solution depends on the stability of the problem. Most practical problems are at least moderately stable, and the true error is then comparable to the error tolerance. To judge the accuracy of the numerical solution, you could reduce tol substantially, e.g., use 0.1 × tol0.1×tol, and solve the problem again. This will usually result in a rather more accurate solution, and the true error of the first integration can be estimated by comparison. Alternatively, a global error assessment can be computed automatically using the parameter errass. Because indirect control of the true error by controlling the local error is generally satisfactory and because both ways of assessing true errors cost twice, or more, the cost of the integration itself, such assessments are used mostly for spot checks, selecting appropriate tolerances for local error control, and exploratory computations.
nag_ode_ivp_rk_range (d02pc) and nag_ode_ivp_rk_onestep (d02pd) each implement three Runge–Kutta formula pairs, and you must select one for the integration. The best choice for method depends on the problem. The order of accuracy is 33, 55 and 88 respectively. As a rule, the smaller tol is, the larger you should take the value of method. If the components thres are small enough that you are effectively specifying relative error control, experience suggests
tol efficient method
10210410-2-10-4 1
10310610-3-10-6 2
10510-5- 3
The overlap in the ranges of tolerances appropriate for a given method merely reflects the dependence of efficiency on the problem being solved. Making tol smaller will normally make the integration more expensive. However, in the range of tolerances appropriate to a method, the increase in cost is modest. There are situations for which one method, or even this kind of code, is a poor choice. You should not specify a very small value for thres(i)thresi, when the iith solution component might vanish. In particular, you should not do this when yi = 0.0yi=0.0. If you do, the code will have to work hard with any value for method to compute significant digits, but method = 1method=1 is a particularly poor choice in this situation. All three methods are inefficient when the problem is ‘stiff’. If it is only mildly stiff, you can solve it with acceptable efficiency with method = 1method=1, but if it is moderately or very stiff, a code designed specifically for such problems will be much more efficient. The higher the order, i.e., the larger the value of method, the more smoothness is required of the solution in order for the method to be efficient.
When assessment of the true (global) error is requested, this error assessment is updated at each step. Its value can be obtained at any time by a call to nag_ode_ivp_rk_errass (d02pz). The code monitors the computation of the global error assessment and reports any doubts it has about the reliability of the results. The assessment scheme requires some smoothness of f(t,y)f(t,y), and it can be deceived if ff is insufficiently smooth. At very crude tolerances the numerical solution can become so inaccurate that it is impossible to continue assessing the accuracy reliably. At very stringent tolerances the effects of finite precision arithmetic can make it impossible to assess the accuracy reliably. The cost of this is roughly twice the cost of the integration itself with method = 2method=2 or 33, and three times with method = 1method=1.
The first step of the integration is critical because it sets the scale of the problem. The integrator will find a starting step size automatically if you set the parameter hstart to 0.00.0. Automatic selection of the first step is so effective that you should normally use it. Nevertheless, you might want to specify a trial value for the first step to be certain that the code recognizes the scale on which phenomena occur near the initial point. Also, automatic computation of the first step size involves some cost, so supplying a good value for this step size will result in a less expensive start. If you are confident that you have a good value, provide it via the parameter hstart.

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:     tstart – double scalar
The initial value of the independent variable, t0t0.
2:     ystart(neq) – double array
neq, the dimension of the array, must satisfy the constraint neq1neq1.
y0y0, the initial values of the solution, yiyi, for i = 1,2,,ni=1,2,,n.
3:     tend – double scalar
The final value of the independent variable, tftf, at which the solution is required. tstart and tend together determine the direction of integration.
Constraint: tendtend must be distinguishable from tstart for the method and the precision of the machine being used.
4:     tol – double scalar
A relative error tolerance.
Constraint: 10.0 × machine precisiontol0.0110.0×machine precisiontol0.01.
5:     thres(neq) – double array
neq, the dimension of the array, must satisfy the constraint neq1neq1.
A vector of thresholds.
Constraint: thres(i)sqrt(σ)thresiσ, where σσ is approximately the smallest possible machine number that can be reciprocated without overflow (see nag_machine_real_safe (x02am)).
6:     method – int64int32nag_int scalar
The Runge–Kutta method to be used.
method = 1method=1
A 2(3)2(3) pair is used.
method = 2method=2
A 4(5)4(5) pair is used.
method = 3method=3
A 7(8)7(8) pair is used.
Constraint: method = 1method=1, 22 or 33.
7:     task – string (length ≥ 1)
Determines whether the usual integration task is to be performed using nag_ode_ivp_rk_range (d02pc) or a more complicated task is to be performed using nag_ode_ivp_rk_onestep (d02pd).
task = 'U'task='U'
nag_ode_ivp_rk_range (d02pc) is to be used for the integration.
task = 'C'task='C'
nag_ode_ivp_rk_onestep (d02pd) is to be used for the integration.
Constraint: task = 'U'task='U' or 'C''C'.
8:     errass – logical scalar
Specifies whether a global error assessment is to be computed with the main integration. errass = trueerrass=true specifies that it is.
Constraint: errass = trueerrass=true or falsefalse.
9:     lenwrk – int64int32nag_int scalar
The dimension of the array work as declared in the (sub)program from which nag_ode_ivp_rk_setup (d02pv) is called. ( lenwrk32 × neqlenwrk32×neq is always sufficient.)
Constraints:

Optional Input Parameters

1:     neq – int64int32nag_int scalar
Default: The dimension of the arrays ystart, thres. (An error is raised if these dimensions are not equal.)
nn, the number of ordinary differential equations in the system to be solved by the integration function.
Constraint: neq1neq1.
2:     hstart – double scalar
A value for the size of the first step in the integration to be attempted. The absolute value of hstart is used with the direction being determined by tstart and tend. The actual first step taken by the integrator may be different to hstart if the underlying algorithm determines that hstart is unsuitable. If hstart = 0.0hstart=0.0 then the size of the first step is computed automatically.
Default: 0.00.0

Input Parameters Omitted from the MATLAB Interface

None.

Output Parameters

1:     work(lenwrk) – double array
Contains information for use by nag_ode_ivp_rk_range (d02pc) or nag_ode_ivp_rk_onestep (d02pd). This must be the same array as supplied to nag_ode_ivp_rk_range (d02pc) or nag_ode_ivp_rk_onestep (d02pd). The contents of this array must remain unchanged between calls.
2:     ifail – int64int32nag_int scalar
ifail = 0ifail=0 unless the function detects an error (see [Error Indicators and Warnings]).

Error Indicators and Warnings

Errors or warnings detected by the function:
  ifail = 1ifail=1
On entry,neq < 1neq<1,
ortend is too close to tstart,
ortol > 0.01tol>0.01 or tol < 10 × machine precisiontol<10×machine precision,
orthres(i) < sqrt(σ)thresi<σ, where σσ is approximately the smallest possible machine number that can be reciprocated without overflow (see nag_machine_real_safe (x02am)),
ormethod1method1, 22 or 33,
ortask'U'task'U' or 'C''C',
orlenwrk is too small.

Accuracy

Not applicable.

Further Comments

If task = 'C'task='C' then the value of the parameter tend may be reset during the integration without the overhead associated with a complete restart; this can be achieved by a call to nag_ode_ivp_rk_reset_tend (d02pw).
It is often the case that a solution component yiyi is of no interest when it is smaller in magnitude than a certain threshold. You can inform the code of this by setting thres(i)thresi to this threshold. In this way you avoid the cost of computing significant digits in yiyi when only the fact that it is smaller than the threshold is of interest. This matter is important when yiyi vanishes, and in particular, when the initial value ystart(i)ystarti vanishes. An appropriate threshold depends on the general size of yiyi in the course of the integration. Physical reasoning may help you select suitable threshold values. If you do not know what to expect of yy, you can find out by a preliminary integration using nag_ode_ivp_rk_range (d02pc) with nominal values of thres. As nag_ode_ivp_rk_range (d02pc) steps from t0t0 towards tftf for each i = 1,2,,ni=1,2,,n it forms ymax(i)ymax(i), the largest magnitude of yiyi computed at any step in the integration so far. Using this you can determine more appropriate values for thres for an accurate integration. You might, for example, take thres(i)thresi to be 10 × machine precision10×machine precision times the final value of ymax(i)ymax(i).

Example

function nag_ode_ivp_rk_setup_example
tstart = 0;
ystart = [0;
     1];
tend = 6.283185307179586;
tol = 0.001;
thres = [1e-08;
     1e-08];
method = int64(1);
task = 'Usual Task';
errass = false;
lenwrk = int64(64);
neq = int64(2);
twant = 0.7853981633974483;
ygot = [0; 0];
ymax = [0; 0];
[work, ifail] = ...
    nag_ode_ivp_rk_setup(tstart, ystart, tend, tol, thres, method, task, errass, lenwrk);
[tgot, ygotOut, ypgot, ymaxOut, workOut, ifail] = ...
    nag_ode_ivp_rk_range(@f, neq, twant, ygot, ymax, work)




function [yp] = f(t, y)
% Evaluate derivative vector.
yp = zeros(2, 1);
yp(1) =  y(2);
yp(2) = -y(1);
 

tgot =

    0.7854


ygotOut =

    0.7069
    0.7068


ypgot =

    0.7058
   -0.7084


ymaxOut =

    0.7069
    1.0000


workOut =

    0.0000
    0.0000
    0.0004
   -0.0005
    0.7049
    0.6924
    0.5988
    0.8006
    0.8006
   -0.5988
    0.8109
    0.5842
    0.5842
   -0.8109
   -0.0373
   -0.0263
   -0.0029
    0.0040
    0.5842
   -0.8109
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0


ifail =

                    0


function d02pv_example
tstart = 0;
ystart = [0;
     1];
tend = 6.283185307179586;
tol = 0.001;
thres = [1e-08;
     1e-08];
method = int64(1);
task = 'Usual Task';
errass = false;
lenwrk = int64(64);
neq = int64(2);
twant = 0.7853981633974483;
ygot = [0; 0];
ymax = [0; 0];
[work, ifail] = ...
    d02pv(tstart, ystart, tend, tol, thres, method, task, errass, lenwrk);
[tgot, ygotOut, ypgot, ymaxOut, workOut, ifail] = ...
    d02pc(@f, neq, twant, ygot, ymax, work)




function [yp] = f(t, y)
% Evaluate derivative vector.
yp = zeros(2, 1);
yp(1) =  y(2);
yp(2) = -y(1);
 

tgot =

    0.7854


ygotOut =

    0.7069
    0.7068


ypgot =

    0.7058
   -0.7084


ymaxOut =

    0.7069
    1.0000


workOut =

    0.0000
    0.0000
    0.0004
   -0.0005
    0.7049
    0.6924
    0.5988
    0.8006
    0.8006
   -0.5988
    0.8109
    0.5842
    0.5842
   -0.8109
   -0.0373
   -0.0263
   -0.0029
    0.0040
    0.5842
   -0.8109
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0


ifail =

                    0



PDF version (NAG web site, 64-bit version, 64-bit version)
Chapter Contents
Chapter Introduction
NAG Toolbox

© The Numerical Algorithms Group Ltd, Oxford, UK. 2009–2013