PUBLIC

RungeKuttaC::RungeKuttaC(void) RungeKuttaC::RungeKuttaC(T &,RealT,RealT,RealT,RealT,BooleanT,int) RungeKuttaC::RungeKuttaC(const RungeKuttaC &) RungeKuttaC::operator=(const RungeKuttaC &) RungeKuttaC::~RungeKuttaC(void) RungeKuttaC::HalveStepsize(RealT &,RealT) RungeKuttaC::DoubleStepsize(BooleanT,RealT) RungeKuttaC::NotClosetoFinalT(RealT,RealT) RungeKuttaC::SetIntstep(RealT) RungeKuttaC::DisplayErrorValues(const T &,const T &) RungeKuttaC::DisplayAnnounceHalving(const RealT &) RungeKuttaC::DisplayAnnounceDoubling(const RealT &) RungeKuttaC::DisplayErrorMessage1(const RealT &) RungeKuttaC::DisplayErrorMessage2(const RealT &) RungeKuttaC::DisplayStep0(const T &,const T &,const T &,const T &) RungeKuttaC::DisplayStep1(const T &,const T &,const T &) RungeKuttaC::DisplayStep2(const T &,const T &,const T &,const T &) RungeKuttaC::DisplayStep3(const T &,const T &,const T &,const T &) RungeKuttaC::DisplayStep4(const T &,const T &,const T &,const T &,const T &) RungeKuttaC::DisplayStep5(const T &,const T &,const T &,const T &,const T &) RungeKuttaC::DisplayInitMessage(const RealT &,const RealT &) RungeKuttaC::DisplaySecondMessage(const RealT &,const RealT &,const RealT &) RungeKuttaC::DisplayValues1(const RealT &,const RealT &,const RealT &) RungeKuttaC::DisplayValues2(const T &,const RealT &,const RealT &,const RealT &)

RungeKuttaC<class T>
A fifth order Runge-Kutta scheme, also called Kutta-Merson scheme.
include	"amma/RungeKutta.hh"
User Level:	Default
Library:	diffeq
Example:	ex_kutta.cc
Section:	Numerical Methods.Differential Equations
In Scope:	std

Comments:
This scheme is for the integration of a set of first order equations, and uses a variable stepsize: starting from a known value (y0), five new values (y1, y2, y3, y4, y5) are calculated using a halve, a third, a sixth and an eight of the current stepsize (h) and the gradients at the calculated values. The difference between the fourth (y4) and the fifth (y5) value is an estimation of the error in the fifth value (y5) compared tot the exact value at a distance h from the old value (y0).

This error in y5 is considered too large if (0.2*Abs(y4-y5) > relerr*Abs(y4)) or (0.2*Abs(y4-y5) > abserr) for all elements of y, in which case the stepsize will be halved. On the other hand, if (0.2*Abs(y4-y5) <= 64*relerr*Abs(y4)) or (0.2*Abs(y4-y5) <= 64*abserr) then the error is much smaller than necessary and the stepsize will be doubled. The initial value for the stepsize ('time step') is given by the user and is equal to the initial value for the stepsize h.

The user can also give a maximum integration time or a number of 'time steps' to be taken. When more than one 'timestep' is to be made, the value for the stepsize h at the end of a 'timestep' is used for the next 'time step'.

----------------------------------------------------------------------
Class RungeKuttaC is made as a template of type T. In order to use this function type T must have:

constuctor, copy constructor, destructor
T operator=(const T& v);
T operator+(const T& v);
T operator-(const T& v);
T operator*(const RealT& d);
void Assign(const T & v);
T Copy(void);
T Abs(void);
BooleanT LargerThan(const T& q, const RealT& threshold)
BooleanT LargerThan(const RealT& val, const RealT& threshold)
BooleanT IsEmpty()
void DisplayVal();

---------------------------------------------------------------------
RungeKuttaC is very similar in form to GNumIntegC, see also the documentation for it. In particular notes on the internal reference.

Class RungeKuttaC is meant to execute one Runge-Kutta integration step. One such step consists of several intermediate steps.
-------------------------------------------------------------------------

RealT step;

initial stepsize, 'time step'

RealT relerr;

maximum relative and absolute error

RealT abserr;

maximum relative and absolute error

BooleanT trace;

trace variable; 0 = no output, 1 = output

T q;

vector with state variables; to be updated

T dqdt;

vector with time derivative of the state variables

int ncuts;

number of stepsize halvings the initial value is equal to the maximum number of stepsize halvings as decided by the user

int neval;

number of function evaluations

RealT h;

current stepsize and its halve, third, sixth and eight part

RealT h2;

current stepsize and its halve, third, sixth and eight part

RealT h3;

current stepsize and its halve, third, sixth and eight part

RealT h6;

current stepsize and its halve, third, sixth and eight part

RealT h8;

current stepsize and its halve, third, sixth and eight part

RungeKuttaC()

RungeKuttaC(T & _q,RealT _step,RealT _h,RealT _relerr,RealT _abserr,BooleanT _tr,int _ncuts)

RungeKuttaC(const RungeKuttaC & qq)

RungeKuttaC & operator=(const RungeKuttaC & qq)

~RungeKuttaC()

Member functions:

-----------------------

int HalveStepsize(RealT & tt,RealT VerySmall)

halves the stepsize

void DoubleStepsize(BooleanT trace,RealT tt)

doubles the stepsize

BooleanT NotClosetoFinalT(RealT tt,RealT FinalT)

checks if the current time is not too close to the final time to double the stepsize

void SetIntstep(RealT step)

sets the integration step

void DisplayErrorValues(const T & err,const T & test)

displays 'err' and 'test'

void DisplayAnnounceHalving(const RealT & time)

void DisplayAnnounceDoubling(const RealT & time)

void DisplayErrorMessage1(const RealT & time)

void DisplayErrorMessage2(const RealT & time)

void DisplayStep0(const T & x0,const T & x1,const T & x2,const T & x3)

void DisplayStep1(const T & x0,const T & x1,const T & x2)

void DisplayStep2(const T & x0,const T & x1,const T & x2,const T & x3)

void DisplayStep3(const T & x0,const T & x1,const T & x2,const T & x3)

void DisplayStep4(const T & x0,const T & x1,const T & x2,const T & x3,const T & x4)

void DisplayStep5(const T & x0,const T & x1,const T & x2,const T & x3,const T & x4)

void DisplayInitMessage(const RealT & x0,const RealT & x1)

void DisplaySecondMessage(const RealT & x0,const RealT & x1,const RealT & x2)

void DisplayValues1(const RealT & x0,const RealT & x1,const RealT & x2)

void DisplayValues2(const T & x0,const RealT & x1,const RealT & x2,const RealT & x3)