Class ODESolver#

Class Documentation#

class ODESolver#

The ODE solver class.

The function of this code is completely explained and documented in the text, Computer solution of ordinary differential equations, the initial value problem by L. F. Shampine and M. K. Gordon. Further details on use of this code are available in Solving ordinary differential equations with ODE, STEP, and INTRP, by L. F. Shampine and M. K. Gordon, SLA-73-1060.

Public Functions

ODESolver(void)#

Create a new ODESolver object with neqn = 0 and no f object.

ODESolver(int n, Derivable *f)#

Create a new solver object with a specific ODE in mind.

Parameters:

  • n – the number of equations

  • f – the derivative function object

int num_eqn(void)#

Get the number of equations this solver was configured for.

Returns:

number of equations

void ode(vector<double> &y, double &t, double tout, double &relerr, double &abserr, int &iflag)#

Integrate a system of neqn first order ODEs.

Abstract

Subroutine ode integrates a system of \(n_{eqn}\) first order ordinary differential equations of the form

\[\frac{d}{dt}y_i = f(t, y_1, y_2, \dots, y_{n})\]

with \(y_i\) given at \(t\).

The subroutine integrates from \(t\) to \(t_{out}\). On return the parameters in the call list are set for continuing the integration. the user has only to define a new value \(t_{out}\) and call ode() again.

The differential equations are actually solved by a suite of codes de(), step1(), and intrp(). ode() allocates virtual storage [as class variables] and calls de() de() is a supervisor which directs the solution. It calls on the routines step1() and intrp() to advance the integration and to interpolate at output points. step1() uses a modified divided difference form of the Adams PECE formulas and local extrapolation. It adjusts the order and step size to control the local error per unit step in a generalized sense. Normally each call to step1() advances the solution one step in the direction of \(t_{out}\). For reasons of efficiency de() integrates beyond \(t_{out}\) internally, though never beyond \(t+10(t_{out}-t)\), and calls intrp() to interpolate the solution at \(t_{out}\). An option is provided to stop the integration at \(t_{out}\) but it should be used only if it is impossible to continue the integration beyond \(t_{out}\).

This code is completely explained and documented in the text, Computer solution of ordinary differential equations, the initial value problem by L. F. Shampine and M. K. Gordon. for usage details, also see SLA-73-1060.

The parameters represent:

  • f – subroutine \(f(t,y,\dot{y})\) to evaluate derivatives \(\dot{y}(i)=dy(i)/dt\)

  • neqn – number of equations to be integrated, \(n_{eqn}\)

  • y[] – solution vector at \(t\)

  • t – independent variable \(t\)

  • tout – point at which solution is desired \(t_{out}\)

  • relerr, abserr – relative and absolute error tolerances for local error test. At each step the code requires \(|\epsilon_{local}| \le |y|\times \epsilon_{rel} + \epsilon_{abs}\) for each component of the local error and solution vectors

  • iflag – indicates status of integration

  • work – [object] to hold information internal to code which is necessary for subsequent calls

First call to ode:

The user must provide storage in his calling program for the arrays in the call list, y(neqn), work(100+21*neqn), iwork(5),declare f in an external statement, supply the subroutine \(f(t,y,\dot{y})\) to evaluate \(dy(i)/dt = \dot{y}(i) = f(t,y(1),y(2),...,y(n_{eqn}))\) and initialize the parameters:

  • neqn – number of equations to be integrated

  • y[] – vector of initial conditions

  • t – starting point of integration

  • tout – point at which solution is desired

  • relerr, abserr – relative and absolute local error tolerances

  • iflag – +1, -1. indicator to initialize the code. Normal input is +1. The user should set iflag==-1 only if it is impossible to continue the integration beyond \(t_{out}\). All parameters except f, neqn and tout may be altered by the code on output so must be variables in the calling program.

Output from ode:

  • neqn – unchanged

  • y[] – solution at \(t\)

  • t – last point reached in integration; normal return has t == tout.

  • tout – unchanged

  • relerr, abserr – normal return has tolerances unchanged. iflag == 3 signals tolerances increased

  • iflag:

    • 2 – normal return. integration reached \(t_{out}\)

    • 3 – integration did not reach \(t_{out}\) because error tolerances too small. relerr, abserr increased appropriately for continuing

    • 4 – integration did not reach \(t_{out}\) because more than 500 steps needed

    • 5 – integration did not reach \(t_{out}\) because equations appear to be stiff

    • 6 – integration did not reach \(t_{out}\) because solution vanished making pure relative error impossible. Must use non-zero abserr to continue.

    • 7 – invalid input parameters (fatal error) the value of iflag is returned negative when the input value is negative and the integration does not reach tout, i.e., -3, -4, -5, -6.

Subsequent calls to ode:

Subroutine returns with all information needed to continue the integration. If the integration reached \(t_{out}\), the user need only define a new \(t_{out}\) and call again. If the integration did not reach \(t_{out}\) and the user wants to continue he just calls again. In the case iflag==6 , the user must also alter the error criterion. The output value of iflag is the appropriate input value for subsequent calls. The only situation in which it should be altered is to stop the integration internally at the new \(t_{out}\), i.e., change output iflag=2 to input iflag=-2. Error tolerances may be changed by the user before continuing. All other parameters must remain unchanged.

Parameters:

  • y – solution vector at t

  • t – independent variable

  • tout – point at which solution is desired

  • relerr – relative error tolerance for local error test

  • abserr – absolute error tolerance for local error test

  • iflag – status of integration

Protected Functions

void de(vector<double> &y, double &t, double tout, double &relerr, double &abserr, int &iflag, bool &phase1, bool &start, int &ns, bool &nornd, int &k, int &kold, int &isnold)#

Differential equation solver.

Abstract

ode merely allocates storage for de to relieve the user of the inconvenience of a long call list and provides the public interface to the DE function.

This code is completely explained and documented in the text, Computer solution of ordinary differential equations, the initial value problem by L. F. Shampine and M. K. Gordon.

Parameters:

  • y – solution vector at T

  • t – independent variable

  • tout – point at which solution is desired

  • relerr – relative error tolerance

  • abserr – absolute error tolerance

  • iflag – status of integration

  • phase1 – work variable

  • start – work variable

  • ns – work variable

  • nornd – work variable

  • k – work variable

  • kold – work variable

  • isnold – work variable

void intrp(double xout, vector<double> &yout, int kold)#

Approximate the solution near x by polymonial.

Abstract

The methods in subroutine step1() approximate the solution near \(x\) by a polynomial. Subroutine intrp() approximates the solution at \(x_{out}\) by evaluating the polynomial there. Information defining this polynomial is passed from step1() so intrp() cannot be used alone.

This code is completely explained and documented in the text, Computer solution of ordinary differential equations, the initial value problem by L. F. Shampine and M. K. Gordon. Further details on use of this code are available in Solving ordinary differential equations with ODE, STEP, and INTRP, by L. F. Shampine and M. K. Gordon, SLA-73-1060.

Input to intrp()

The user provides storage in the calling program for the arrays in the call list - dimension y(neqn),yout(neqn),ypout(neqn),phi(neqn,16),psi(12)

and defines - xout – point at which solution is desired.

The remaining parameters are defined by a previous call to step1() and should be passed to intrp() unchanged.

Output from intrp()

  • yout – solution at \(x_{out}\)

  • ypout – derivative of solution at \(x_{out}\)

the remaining parameters are returned unaltered from their input values. Integration with step1() may be continued.

Parameters:

  • xout – points at which solution is known

  • yout – solution at xout

  • kold – pass from sansmic::step1 - unchanged

void step1(double &eps, bool &start, int &k, int &kold, bool &crash, bool &phase1, int &ns, bool &nornd)#

Take a single solver step, used indirectly through ODE.

Abstract

Subroutine step1() is normally used indirectly through subroutine ode(). Because ode() suffices for most problems and is much easier to use, using it should be considered before using step1() alone.

Subroutine step1() integrates a system of \(n_{eqn}\) first order ordinary differential equations one step, normally from \(x\) to \(x+h\), using a modified divided difference form of the Adams PECE formulas. Local extrapolation is used to improve absolute stability and accuracy. The code adjusts its order and step size to control the local error per unit step in a generalized sense. special devices are included to control roundoff error and to detect when the user is requesting too much accuracy.

This code is completely explained and documented in the text, Computer solution of ordinary differential equations, the initial value problem by L. F. Shampine and M. K. Gordon. Further details on use of this code are available in Solving ordinary differential equations with ODE, STEP, and INTRP, by L. F. Shampine and M. K. Gordon, SLA-73-1060.

The parameters represent: - f – subroutine to evaluate derivatives - neqn – number of equations to be integrated - y[] – solution vector at \(x\) - x – independent variable - h – appropriate step size for next step. normally determined by code - eps – local error tolerance - wt[] – vector of weights for error criterion - start – logical variable set true for first step, false otherwise - hold – step size used for last successful step - k – appropriate order for next step (determined by code) - kold – order used for last successful step - crash – logical variable set true when no step can be taken, else false - yp[] – derivative of solution vector at \(x\) after successful step

The arrays phi, psi are required for the interpolation subroutine intrp(). The array p is internal to the code. The remaining nine variables and arrays are included in the call list only to eliminate local retention of variables between calls.

Input to step1()
First call

The user must provide storage in his calling program for all arrays in the call list, namely

dimension
y(neqn),wt(neqn),phi(neqn,16),p(neqn),yp(neqn),psi(12), 1
alpha(12),beta(12),sig(13),v(12),w(12),g(13)

The user must also declare start, crash, phase1 and nornd logical variables and f an external subroutine, supply the subroutine f(x,y,yp) to evaluate

\[\frac{d}{dx} y_i = \dot{y}_i = f(x, y_1, y_2, \dots, y_n)\]

and initialize only the following parameters. (Note: this is provided by the PlumeRise class for the specific SANSMIC case.)

  • neqn – number of equations to be integrated

  • y[] – vector of initial values of dependent variables

  • x – initial value of the independent variable

  • h – nominal step size indicating direction of integration and maximum size of step. must be variable

  • eps – local error tolerance per step. must be variable

  • wt[] – vector of non-zero weights for error criterion

  • start – true

step1 requires that the l2 norm of the vector with components local error(l)/wt(l) be less than eps for a successful step. the array wt allows the user to specify an error test appropriate for their problem. For example,

wt(l) = 1.0

specifies absolute error,

wt(l) = abs(y(l))

error relative to the most recent value of the l-th component of the solution, wt(l) = abs(yp(l)) error relative to the most recent value of the l-th component of the derivative, wt(l) = amax1(wt(l),abs(y(l))) error relative to the largest magnitude of l-th component obtained so far, wt(l) abs(y(l))*relerr/eps + abserr/eps specifies a mixed relative-absolute test where relerr is relative error, abserr is absolute error and eps = amax1(relerr,abserr) .

Subsequent calls

Subroutine step1() is designed so that all information needed to continue the integration, including the step size h and the order k , is returned with each step. with the exception of the step size, the error tolerance, and the weights, none of the parameters should be altered. the array wt must be updated after each step to maintain relative error tests like those above. normally the integration is continued just beyond the desired endpoint and the solution interpolated there with subroutine intrp . if it is impossible to integrate beyond the endpoint, the step size may be reduced to hit the endpoint since the code will not take a step larger than the h input. changing the direction of integration, i.e., the sign of h , requires the user set start = .true. before calling step1 again. this is the only situation in which start should be altered.

Output from STEP1
Successful step

the subroutine returns after each successful step with start and crash set .false. . x represents the independent variable advanced one step of length hold from its value on input and y the solution vector at the new value of x . all other parameters represent information corresponding to the new x needed to continue the integration.

Unsuccessful step

when the error tolerance is too small for the machine precision, the subroutine returns without taking a step and crash = .true. . an appropriate step size and error tolerance for continuing are estimated and all other information is restored as upon input before returning. to continue with the larger tolerance, the user just calls the code again. a restart is neither required nor desirable.

Parameters:

  • eps – local error tolerance

  • start – logical variable set true for first step

  • k – appropriate order for next step

  • kold – order used for last successful step

  • crash – logical variable set true when no step can be taken

  • phase1 – elimainate local retention of variables

  • ns – elimainate local retention of variables

  • nornd – elimainate local retention of variables

Protected Attributes

int neqn#

the number of equations

Derivable *dy_dx#

the derivative function object