C++ version of DUMKA3 is a template function dumka3.cpp.
Idea to use templates was taken from IML++ library. In order to use
templates you need to have classes VectorFunction and Vector
-
VectorFunction - implements VectorFunctionBase - interface with
three virtual functions, which have to be implemented
-
void compute(Vector& y,
double t, Vector& result) - calculates right hand side of the ODE
y'=f(y,t)
-
double cour(Vector& y,
double t) - used if "isComputeEigenValue=false" to
-
calculate
h_{euler} - step which you would use in explicit Euler method.
-
If
"isComputeEigenValue=true" internal eigenvalue
function will be used to evaluate h_{euler}=2/max(eigenvalue).
-
Element GetNonZeroElement()
- can be used if "isComputeEigenValue=true" in eigenvalue
function, if
-
iterations
of power method have to be restarted, but y,f(y) equals to zero,
-
in this
case program need to initialize random non-zero vector, because
-
vector's
Element is unknown in advance
-
for
example if Vector contains "double"(s) you can return
rand()%100/100.
-
but if
Vector contains "Point"(s) you can
-
return
Point(rand()%100/100.,rand()%100/100.,rand()%100/100.)
-
-
template <class Vector,
class Element>
-
class VectorFunctionBase {
-
public:
-
virtual void
compute(Vector& y, double t, Vector& result) = 0;
-
virtual double
cour(Vector& y, double t) = 0;
-
virtual Element
GetNonZeroElement() = 0;
-
};
-
Vector class should have function size() and operator[]
(for example class vector<double> has these functions).
-
Each Element of the Vector should have operator* (for
example if your Vector has elements "double" everything is
OK, but if you use "Point" you should define operator*)
This is the first C++ version of dumka3, and as everything new it will
be changed and modified in a future, and your
suggestions/complains/questions would be very welcome.
Both fortran version of dumka3.for and C++
version of dumka3.cpp require additional
function which determine maximal step of explicit Euler method. It can
be evaluated by value cou=2/M, where M is upper boundary of maximal
eigen value. For example if you solve ODE y'=-5*y, then take
cou=2/5 or something smaller, for example 2/5.1. You can use also norm
of jacobian of the right hand side function.