849 lines
32 KiB
MQL5
849 lines
32 KiB
MQL5
//+------------------------------------------------------------------+
|
|
//| diffequations.mqh |
|
|
//| Copyright 2003-2022 Sergey Bochkanov (ALGLIB project) |
|
|
//| Copyright 2012-2023, MetaQuotes Ltd. |
|
|
//| https://www.mql5.com |
|
|
//+------------------------------------------------------------------+
|
|
//| Implementation of ALGLIB library in MetaQuotes Language 5 |
|
|
//| |
|
|
//| The features of the library include: |
|
|
//| - Linear algebra (direct algorithms, EVD, SVD) |
|
|
//| - Solving systems of linear and non-linear equations |
|
|
//| - Interpolation |
|
|
//| - Optimization |
|
|
//| - FFT (Fast Fourier Transform) |
|
|
//| - Numerical integration |
|
|
//| - Linear and nonlinear least-squares fitting |
|
|
//| - Ordinary differential equations |
|
|
//| - Computation of special functions |
|
|
//| - Descriptive statistics and hypothesis testing |
|
|
//| - Data analysis - classification, regression |
|
|
//| - Implementing linear algebra algorithms, interpolation, etc. |
|
|
//| in high-precision arithmetic (using MPFR) |
|
|
//| |
|
|
//| This file is free software; you can redistribute it and/or |
|
|
//| modify it under the terms of the GNU General Public License as |
|
|
//| published by the Free Software Foundation (www.fsf.org); either |
|
|
//| version 2 of the License, or (at your option) any later version. |
|
|
//| |
|
|
//| This program is distributed in the hope that it will be useful, |
|
|
//| but WITHOUT ANY WARRANTY; without even the implied warranty of |
|
|
//| MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the |
|
|
//| GNU General Public License for more details. |
|
|
//+------------------------------------------------------------------+
|
|
#include "matrix.mqh"
|
|
#include "alglibinternal.mqh"
|
|
//+------------------------------------------------------------------+
|
|
//| Auxiliary class for CODESolver |
|
|
//+------------------------------------------------------------------+
|
|
class CODESolverState
|
|
{
|
|
public:
|
|
int m_n;
|
|
int m_m;
|
|
double m_xscale;
|
|
double m_h;
|
|
double m_eps;
|
|
bool m_fraceps;
|
|
int m_repterminationtype;
|
|
int m_repnfev;
|
|
int m_solvertype;
|
|
bool m_needdy;
|
|
double m_x;
|
|
RCommState m_rstate;
|
|
//--- arrays
|
|
double m_yc[];
|
|
double m_escale[];
|
|
double m_xg[];
|
|
double m_y[];
|
|
double m_dy[];
|
|
double m_yn[];
|
|
double m_yns[];
|
|
double m_rka[];
|
|
double m_rkc[];
|
|
double m_rkcs[];
|
|
//--- matrices
|
|
CMatrixDouble m_ytbl;
|
|
CMatrixDouble m_rkb;
|
|
CMatrixDouble m_rkk;
|
|
|
|
public:
|
|
CODESolverState(void);
|
|
~CODESolverState(void) {}
|
|
|
|
void Copy(CODESolverState &obj);
|
|
};
|
|
//+------------------------------------------------------------------+
|
|
//| Constructor |
|
|
//+------------------------------------------------------------------+
|
|
CODESolverState::CODESolverState(void)
|
|
{
|
|
m_n=0;
|
|
m_m=0;
|
|
m_xscale=0;
|
|
m_h=0;
|
|
m_eps=0;
|
|
m_fraceps=false;
|
|
m_repterminationtype=0;
|
|
m_repnfev=0;
|
|
m_solvertype=0;
|
|
m_needdy=false;
|
|
m_x=0;
|
|
}
|
|
//+------------------------------------------------------------------+
|
|
//| Copy |
|
|
//+------------------------------------------------------------------+
|
|
void CODESolverState::Copy(CODESolverState &obj)
|
|
{
|
|
//--- copy variables
|
|
m_n=obj.m_n;
|
|
m_m=obj.m_m;
|
|
m_xscale=obj.m_xscale;
|
|
m_h=obj.m_h;
|
|
m_eps=obj.m_eps;
|
|
m_fraceps=obj.m_fraceps;
|
|
m_repterminationtype=obj.m_repterminationtype;
|
|
m_repnfev=obj.m_repnfev;
|
|
m_solvertype=obj.m_solvertype;
|
|
m_needdy=obj.m_needdy;
|
|
m_x=obj.m_x;
|
|
m_rstate.Copy(obj.m_rstate);
|
|
//--- copy arrays
|
|
ArrayCopy(m_yc,obj.m_yc);
|
|
ArrayCopy(m_escale,obj.m_escale);
|
|
ArrayCopy(m_xg,obj.m_xg);
|
|
ArrayCopy(m_y,obj.m_y);
|
|
ArrayCopy(m_dy,obj.m_dy);
|
|
ArrayCopy(m_yn,obj.m_yn);
|
|
ArrayCopy(m_yns,obj.m_yns);
|
|
ArrayCopy(m_rka,obj.m_rka);
|
|
ArrayCopy(m_rkc,obj.m_rkc);
|
|
ArrayCopy(m_rkcs,obj.m_rkcs);
|
|
//--- copy matrices
|
|
m_ytbl=obj.m_ytbl;
|
|
m_rkb=obj.m_rkb;
|
|
m_rkk=obj.m_rkk;
|
|
}
|
|
//+------------------------------------------------------------------+
|
|
//| This class is a shell for class CODESolverState |
|
|
//+------------------------------------------------------------------+
|
|
class CODESolverStateShell
|
|
{
|
|
private:
|
|
CODESolverState m_innerobj;
|
|
|
|
public:
|
|
//--- constructors, destructor
|
|
CODESolverStateShell(void) {}
|
|
CODESolverStateShell(CODESolverState &obj) { m_innerobj.Copy(obj); }
|
|
~CODESolverStateShell(void) {}
|
|
//--- methods
|
|
bool GetNeedDY(void);
|
|
void SetNeedDY(const bool b);
|
|
double GetX(void);
|
|
void SetX(const double d);
|
|
CODESolverState *GetInnerObj(void);
|
|
};
|
|
//+------------------------------------------------------------------+
|
|
//| Returns the value of the variable needdy |
|
|
//+------------------------------------------------------------------+
|
|
bool CODESolverStateShell::GetNeedDY(void)
|
|
{
|
|
return(m_innerobj.m_needdy);
|
|
}
|
|
//+------------------------------------------------------------------+
|
|
//| Changing the value of the variable needdy |
|
|
//+------------------------------------------------------------------+
|
|
void CODESolverStateShell::SetNeedDY(const bool b)
|
|
{
|
|
m_innerobj.m_needdy=b;
|
|
}
|
|
//+------------------------------------------------------------------+
|
|
//| Returns the value of the variable x |
|
|
//+------------------------------------------------------------------+
|
|
double CODESolverStateShell::GetX(void)
|
|
{
|
|
return(m_innerobj.m_x);
|
|
}
|
|
//+------------------------------------------------------------------+
|
|
//| Changing the value of the variable x |
|
|
//+------------------------------------------------------------------+
|
|
void CODESolverStateShell::SetX(const double d)
|
|
{
|
|
m_innerobj.m_x=d;
|
|
}
|
|
//+------------------------------------------------------------------+
|
|
//| Return object of class |
|
|
//+------------------------------------------------------------------+
|
|
CODESolverState *CODESolverStateShell::GetInnerObj(void)
|
|
{
|
|
return(GetPointer(m_innerobj));
|
|
}
|
|
//+------------------------------------------------------------------+
|
|
//| Auxiliary class for CODESolver |
|
|
//+------------------------------------------------------------------+
|
|
class CODESolverReport
|
|
{
|
|
public:
|
|
//--- class variables
|
|
int m_nfev;
|
|
int m_terminationtype;
|
|
//--- constructor, destructor
|
|
CODESolverReport(void) { ZeroMemory(this); }
|
|
~CODESolverReport(void) {}
|
|
//--- copy
|
|
void Copy(CODESolverReport &obj);
|
|
};
|
|
//+------------------------------------------------------------------+
|
|
//| Copy |
|
|
//+------------------------------------------------------------------+
|
|
void CODESolverReport::Copy(CODESolverReport &obj)
|
|
{
|
|
//--- copy variables
|
|
m_nfev=obj.m_nfev;
|
|
m_terminationtype=obj.m_terminationtype;
|
|
}
|
|
//+------------------------------------------------------------------+
|
|
//| This class is a shell for class CODESolverReport |
|
|
//+------------------------------------------------------------------+
|
|
class CODESolverReportShell
|
|
{
|
|
private:
|
|
CODESolverReport m_innerobj;
|
|
|
|
public:
|
|
//--- constructor, destructor
|
|
CODESolverReportShell(void) {}
|
|
CODESolverReportShell(CODESolverReport &obj) { m_innerobj.Copy(obj); }
|
|
~CODESolverReportShell(void) {}
|
|
//--- methods
|
|
int GetNFev(void);
|
|
void SetNFev(const int i);
|
|
int GetTerminationType(void);
|
|
void SetTerminationType(const int i);
|
|
CODESolverReport *GetInnerObj(void);
|
|
};
|
|
//+------------------------------------------------------------------+
|
|
//| Returns the value of the variable nfev |
|
|
//+------------------------------------------------------------------+
|
|
int CODESolverReportShell::GetNFev(void)
|
|
{
|
|
return(m_innerobj.m_nfev);
|
|
}
|
|
//+------------------------------------------------------------------+
|
|
//| Changing the value of the variable nfev |
|
|
//+------------------------------------------------------------------+
|
|
void CODESolverReportShell::SetNFev(const int i)
|
|
{
|
|
m_innerobj.m_nfev=i;
|
|
}
|
|
//+------------------------------------------------------------------+
|
|
//| Returns the value of the variable terminationtype |
|
|
//+------------------------------------------------------------------+
|
|
int CODESolverReportShell::GetTerminationType(void)
|
|
{
|
|
return(m_innerobj.m_terminationtype);
|
|
}
|
|
//+------------------------------------------------------------------+
|
|
//| Changing the value of the variable terminationtype |
|
|
//+------------------------------------------------------------------+
|
|
void CODESolverReportShell::SetTerminationType(const int i)
|
|
{
|
|
m_innerobj.m_terminationtype=i;
|
|
}
|
|
//+------------------------------------------------------------------+
|
|
//| Return object of class |
|
|
//+------------------------------------------------------------------+
|
|
CODESolverReport *CODESolverReportShell::GetInnerObj(void)
|
|
{
|
|
return(GetPointer(m_innerobj));
|
|
}
|
|
//+------------------------------------------------------------------+
|
|
//| Solution of ordinary differential equations |
|
|
//+------------------------------------------------------------------+
|
|
class CODESolver
|
|
{
|
|
public:
|
|
//--- class constants
|
|
static const double m_odesolvermaxgrow;
|
|
static const double m_odesolvermaxshrink;
|
|
//--- public methods
|
|
static void ODESolverRKCK(double &y[],const int n,double &x[],const int m,const double eps,const double h,CODESolverState &state);
|
|
static void ODESolverResults(CODESolverState &state,int &m,double &xtbl[],CMatrixDouble &ytbl,CODESolverReport &rep);
|
|
static bool ODESolverIteration(CODESolverState &state);
|
|
|
|
private:
|
|
//--- private method
|
|
static void ODESolverInit(int solvertype,double &y[],const int n,double &x[],const int m,const double eps,double h,CODESolverState &state);
|
|
//--- auxiliary functions for ODESolverIteration
|
|
static void Func_lbl_rcomm(CODESolverState &state,int n,int m,int i,int j,int k,int klimit,bool gridpoint,double xc,double v,double h,double h2,double err,double maxgrowpow);
|
|
static bool Func_lbl_6(CODESolverState &state,int &n,int &m,int &i,int &j,int &k,int &klimit,bool &gridpoint,double &xc,double &v,double &h,double &h2,double &err,double &maxgrowpow);
|
|
static bool Func_lbl_8(CODESolverState &state,int &n,int &m,int &i,int &j,int &k,int &klimit,bool &gridpoint,double &xc,double &v,double &h,double &h2,double &err,double &maxgrowpow);
|
|
static bool Func_lbl_10(CODESolverState &state,int &n,int &m,int &i,int &j,int &k,int &klimit,bool &gridpoint,double &xc,double &v,double &h,double &h2,double &err,double &maxgrowpow);
|
|
};
|
|
//+------------------------------------------------------------------+
|
|
//| Initialize constants |
|
|
//+------------------------------------------------------------------+
|
|
const double CODESolver::m_odesolvermaxgrow=3.0;
|
|
const double CODESolver::m_odesolvermaxshrink=10.0;
|
|
//+------------------------------------------------------------------+
|
|
//| Cash-Karp adaptive ODE solver. |
|
|
//| This subroutine solves ODE Y'=f(Y,x) with initial conditions |
|
|
//| Y(xs)=Ys (here Y may be single variable or vector of N variables)|
|
|
//| INPUT PARAMETERS: |
|
|
//| Y - initial conditions, array[0..N-1]. |
|
|
//| contains values of Y[] at X[0] |
|
|
//| N - system size |
|
|
//| X - points at which Y should be tabulated, |
|
|
//| array[0..M-1] integrations starts at X[0], ends |
|
|
//| at X[M-1], intermediate values at X[i] are |
|
|
//| returned too. |
|
|
//| SHOULD BE ORDERED BY ASCENDING OR BY DESCENDING!!|
|
|
//| M - number of intermediate points + first point + |
|
|
//| last point: |
|
|
//| * M>2 means that you need both Y(X[M-1]) and M-2 |
|
|
//| values at intermediate points |
|
|
//| * M=2 means that you want just to integrate from |
|
|
//| X[0] to X[1] and don't interested in |
|
|
//| intermediate values. |
|
|
//| * M=1 means that you don't want to integrate :) |
|
|
//| it is degenerate case, but it will be handled |
|
|
//| correctly. |
|
|
//| * M<1 means error |
|
|
//| Eps - tolerance (absolute/relative error on each step |
|
|
//| will be less than Eps). When passing: |
|
|
//| * Eps>0, it means desired ABSOLUTE error |
|
|
//| * Eps<0, it means desired RELATIVE error. |
|
|
//| Relative errors are calculated with respect to |
|
|
//| maximum values of Y seen so far. Be careful to |
|
|
//| use this criterion when starting from Y[] that |
|
|
//| are close to zero. |
|
|
//| H - initial step lenth, it will be adjusted |
|
|
//| automatically after the first step. If H=0, step |
|
|
//| will be selected automatically (usualy it will |
|
|
//| be equal to 0.001 of min(x[i]-x[j])). |
|
|
//| OUTPUT PARAMETERS |
|
|
//| State - structure which stores algorithm state between |
|
|
//| subsequent calls of OdeSolverIteration. Used |
|
|
//| for reverse communication. This structure should |
|
|
//| be passed to the OdeSolverIteration subroutine. |
|
|
//| SEE ALSO |
|
|
//| AutoGKSmoothW, AutoGKSingular, AutoGKIteration, AutoGKResults|
|
|
//+------------------------------------------------------------------+
|
|
void CODESolver::ODESolverRKCK(double &y[],const int n,double &x[],
|
|
const int m,const double eps,
|
|
const double h,CODESolverState &state)
|
|
{
|
|
//--- check
|
|
if(!CAp::Assert(n>=1,"ODESolverRKCK: N<1!"))
|
|
return;
|
|
//--- check
|
|
if(!CAp::Assert(m>=1,"ODESolverRKCK: M<1!"))
|
|
return;
|
|
//--- check
|
|
if(!CAp::Assert(CAp::Len(y)>=n,"ODESolverRKCK: Length(Y)<N!"))
|
|
return;
|
|
//--- check
|
|
if(!CAp::Assert(CAp::Len(x)>=m,"ODESolverRKCK: Length(X)<M!"))
|
|
return;
|
|
//--- check
|
|
if(!CAp::Assert(CApServ::IsFiniteVector(y,n),"ODESolverRKCK: Y contains infinite or NaN values!"))
|
|
return;
|
|
//--- check
|
|
if(!CAp::Assert(CApServ::IsFiniteVector(x,m),"ODESolverRKCK: Y contains infinite or NaN values!"))
|
|
return;
|
|
//--- check
|
|
if(!CAp::Assert(CMath::IsFinite(eps),"ODESolverRKCK: Eps is not finite!"))
|
|
return;
|
|
//--- check
|
|
if(!CAp::Assert(eps!=0.0,"ODESolverRKCK: Eps is zero!"))
|
|
return;
|
|
//--- check
|
|
if(!CAp::Assert(CMath::IsFinite(h),"ODESolverRKCK: H is not finite!"))
|
|
return;
|
|
//--- function call
|
|
ODESolverInit(0,y,n,x,m,eps,h,state);
|
|
}
|
|
//+------------------------------------------------------------------+
|
|
//| ODE solver results |
|
|
//| Called after OdeSolverIteration returned False. |
|
|
//| INPUT PARAMETERS: |
|
|
//| State - algorithm state (used by OdeSolverIteration). |
|
|
//| OUTPUT PARAMETERS: |
|
|
//| M - number of tabulated values, M>=1 |
|
|
//| XTbl - array[0..M-1], values of X |
|
|
//| YTbl - array[0..M-1,0..N-1], values of Y in X[i] |
|
|
//| Rep - solver report: |
|
|
//| * Rep.TerminationType completetion code: |
|
|
//| * -2 X is not ordered by |
|
|
//| ascending/descending or there are |
|
|
//| non-distinct X[], i.e. X[i]=X[i+1] |
|
|
//| * -1 incorrect parameters were specified |
|
|
//| * 1 task has been solved |
|
|
//| * Rep.NFEV contains number of function |
|
|
//| calculations |
|
|
//+------------------------------------------------------------------+
|
|
void CODESolver::ODESolverResults(CODESolverState &state,int &m,
|
|
double &xtbl[],CMatrixDouble &ytbl,
|
|
CODESolverReport &rep)
|
|
{
|
|
//--- create variables
|
|
double v=0;
|
|
int i=0;
|
|
int i_=0;
|
|
//--- initialization
|
|
m=0;
|
|
rep.m_terminationtype=state.m_repterminationtype;
|
|
//--- check
|
|
if(rep.m_terminationtype>0)
|
|
{
|
|
//--- change values
|
|
m=state.m_m;
|
|
rep.m_nfev=state.m_repnfev;
|
|
//--- allocation
|
|
ArrayResize(xtbl,state.m_m);
|
|
v=state.m_xscale;
|
|
//--- calculation
|
|
for(i_=0; i_<=state.m_m-1; i_++)
|
|
xtbl[i_]=v*state.m_xg[i_];
|
|
//--- allocation
|
|
ytbl.Resize(state.m_m,state.m_n);
|
|
for(i=0; i<=state.m_m-1; i++)
|
|
{
|
|
for(i_=0; i_<=state.m_n-1; i_++)
|
|
ytbl.Set(i,i_,state.m_ytbl[i][i_]);
|
|
}
|
|
}
|
|
else
|
|
rep.m_nfev=0;
|
|
}
|
|
//+------------------------------------------------------------------+
|
|
//| Internal initialization subroutine |
|
|
//+------------------------------------------------------------------+
|
|
void CODESolver::ODESolverInit(int solvertype,double &y[],const int n,
|
|
double &x[],const int m,const double eps,
|
|
double h,CODESolverState &state)
|
|
{
|
|
//--- create variables
|
|
int i=0;
|
|
double v=0;
|
|
int i_=0;
|
|
//--- Prepare RComm
|
|
state.m_rstate.ia.Resize(6);
|
|
ArrayResizeAL(state.m_rstate.ba,1);
|
|
state.m_rstate.ra.Resize(6);
|
|
state.m_rstate.stage=-1;
|
|
state.m_needdy=false;
|
|
//--- check parameters.
|
|
if((n<=0 || m<1) || eps==0.0)
|
|
{
|
|
state.m_repterminationtype=-1;
|
|
return;
|
|
}
|
|
//--- check
|
|
if(h<0.0)
|
|
h=-h;
|
|
//--- quick exit if necessary.
|
|
//--- after this block we assume that M>1
|
|
if(m==1)
|
|
{
|
|
//--- change values
|
|
state.m_repnfev=0;
|
|
state.m_repterminationtype=1;
|
|
state.m_ytbl.Resize(1,n);
|
|
for(i_=0; i_<=n-1; i_++)
|
|
state.m_ytbl.Set(0,i_,y[i_]);
|
|
//--- allocation
|
|
ArrayResize(state.m_xg,m);
|
|
for(i_=0; i_<=m-1; i_++)
|
|
state.m_xg[i_]=x[i_];
|
|
//--- exit the function
|
|
return;
|
|
}
|
|
//--- check again: correct order of X[]
|
|
if(x[1]==x[0])
|
|
{
|
|
state.m_repterminationtype=-2;
|
|
return;
|
|
}
|
|
for(i=1; i<=m-1; i++)
|
|
{
|
|
//--- check
|
|
if((x[1]>x[0] && x[i]<=x[i-1]) || (x[1]<x[0] && x[i]>=x[i-1]))
|
|
{
|
|
state.m_repterminationtype=-2;
|
|
return;
|
|
}
|
|
}
|
|
//--- auto-select H if necessary
|
|
if(h==0.0)
|
|
{
|
|
v=MathAbs(x[1]-x[0]);
|
|
for(i=2; i<=m-1; i++)
|
|
v=MathMin(v,MathAbs(x[i]-x[i-1]));
|
|
h=0.001*v;
|
|
}
|
|
//--- store parameters
|
|
state.m_n=n;
|
|
state.m_m=m;
|
|
state.m_h=h;
|
|
state.m_eps=MathAbs(eps);
|
|
state.m_fraceps=eps<0.0;
|
|
//--- allocation
|
|
ArrayResize(state.m_xg,m);
|
|
for(i_=0; i_<=m-1; i_++)
|
|
state.m_xg[i_]=x[i_];
|
|
//--- check
|
|
if(x[1]>x[0])
|
|
state.m_xscale=1;
|
|
else
|
|
{
|
|
state.m_xscale=-1;
|
|
for(i_=0; i_<=m-1; i_++)
|
|
state.m_xg[i_]=-1*state.m_xg[i_];
|
|
}
|
|
//--- allocation
|
|
ArrayResize(state.m_yc,n);
|
|
for(i_=0; i_<=n-1; i_++)
|
|
state.m_yc[i_]=y[i_];
|
|
//--- change values
|
|
state.m_solvertype=solvertype;
|
|
state.m_repterminationtype=0;
|
|
//--- Allocate arrays
|
|
ArrayResize(state.m_y,n);
|
|
ArrayResize(state.m_dy,n);
|
|
}
|
|
//+------------------------------------------------------------------+
|
|
//| Iterative method |
|
|
//+------------------------------------------------------------------+
|
|
bool CODESolver::ODESolverIteration(CODESolverState &state)
|
|
{
|
|
//--- create variables
|
|
int n=0;
|
|
int m=0;
|
|
int i=0;
|
|
int j=0;
|
|
int k=0;
|
|
double xc=0;
|
|
double v=0;
|
|
double h=0;
|
|
double h2=0;
|
|
bool gridpoint;
|
|
double err=0;
|
|
double maxgrowpow=0;
|
|
int klimit=0;
|
|
int i_=0;
|
|
//--- This code initializes locals by:
|
|
//--- * random values determined during code
|
|
//--- generation - on first subroutine call
|
|
//--- * values from previous call - on subsequent calls
|
|
if(state.m_rstate.stage>=0)
|
|
{
|
|
//--- initialization
|
|
n=state.m_rstate.ia[0];
|
|
m=state.m_rstate.ia[1];
|
|
i=state.m_rstate.ia[2];
|
|
j=state.m_rstate.ia[3];
|
|
k=state.m_rstate.ia[4];
|
|
klimit=state.m_rstate.ia[5];
|
|
gridpoint=state.m_rstate.ba[0];
|
|
xc=state.m_rstate.ra[0];
|
|
v=state.m_rstate.ra[1];
|
|
h=state.m_rstate.ra[2];
|
|
h2=state.m_rstate.ra[3];
|
|
err=state.m_rstate.ra[4];
|
|
maxgrowpow=state.m_rstate.ra[5];
|
|
}
|
|
else
|
|
{
|
|
//--- initialization
|
|
n=-983;
|
|
m=-989;
|
|
i=-834;
|
|
j=900;
|
|
k=-287;
|
|
klimit=364;
|
|
gridpoint=false;
|
|
xc=-338;
|
|
v=-686;
|
|
h=912;
|
|
h2=585;
|
|
err=497;
|
|
maxgrowpow=-271;
|
|
}
|
|
//--- check
|
|
if(state.m_rstate.stage==0)
|
|
{
|
|
//--- change values
|
|
state.m_needdy=false;
|
|
state.m_repnfev=state.m_repnfev+1;
|
|
v=h*state.m_xscale;
|
|
for(i_=0; i_<=n-1; i_++)
|
|
state.m_rkk.Set(k,i_,v*state.m_dy[i_]);
|
|
//--- update YN/YNS
|
|
v=state.m_rkc[k];
|
|
for(i_=0; i_<=n-1; i_++)
|
|
state.m_yn[i_]=state.m_yn[i_]+v*state.m_rkk[k][i_];
|
|
v=state.m_rkcs[k];
|
|
for(i_=0; i_<=n-1; i_++)
|
|
state.m_yns[i_]=state.m_yns[i_]+v*state.m_rkk[k][i_];
|
|
k=k+1;
|
|
return(Func_lbl_8(state,n,m,i,j,k,klimit,gridpoint,xc,v,h,h2,err,maxgrowpow));
|
|
}
|
|
//--- Routine body
|
|
//--- prepare
|
|
if(state.m_repterminationtype!=0)
|
|
return(false);
|
|
//--- change values
|
|
n=state.m_n;
|
|
m=state.m_m;
|
|
h=state.m_h;
|
|
maxgrowpow=MathPow(m_odesolvermaxgrow,5);
|
|
state.m_repnfev=0;
|
|
//--- some preliminary checks for internal errors
|
|
//--- after this we assume that H>0 and M>1
|
|
if(!CAp::Assert(state.m_h>0.0,"ODESolver: internal error"))
|
|
return(false);
|
|
//--- check
|
|
if(!CAp::Assert(m>1,"ODESolverIteration: internal error"))
|
|
return(false);
|
|
//--- choose solver
|
|
if(state.m_solvertype!=0)
|
|
return(false);
|
|
//--- Cask-Karp solver
|
|
//--- Prepare coefficients table.
|
|
//--- Check it for errors
|
|
ArrayResize(state.m_rka,6);
|
|
//--- calculation
|
|
state.m_rka[0]=0;
|
|
state.m_rka[1]=1.0/5.0;
|
|
state.m_rka[2]=3.0/10.0;
|
|
state.m_rka[3]=3.0/5.0;
|
|
state.m_rka[4]=1;
|
|
state.m_rka[5]=7.0/8.0;
|
|
state.m_rkb.Resize(6,5);
|
|
state.m_rkb.Set(1,0,1.0/5.0);
|
|
state.m_rkb.Set(2,0,3.0/40.0);
|
|
state.m_rkb.Set(2,1,9.0/40.0);
|
|
state.m_rkb.Set(3,0,3.0/10.0);
|
|
state.m_rkb.Set(3,1,-(9.0/10.0));
|
|
state.m_rkb.Set(3,2,6.0/5.0);
|
|
state.m_rkb.Set(4,0,-(11.0/54.0));
|
|
state.m_rkb.Set(4,1,5.0/2.0);
|
|
state.m_rkb.Set(4,2,-(70.0/27.0));
|
|
state.m_rkb.Set(4,3,35.0/27.0);
|
|
state.m_rkb.Set(5,0,1631.0/55296.0);
|
|
state.m_rkb.Set(5,1,175.0/512.0);
|
|
state.m_rkb.Set(5,2,575.0/13824.0);
|
|
state.m_rkb.Set(5,3,44275.0/110592.0);
|
|
state.m_rkb.Set(5,4,253.0/4096.0);
|
|
//--- allocation
|
|
ArrayResize(state.m_rkc,6);
|
|
//--- calculation
|
|
state.m_rkc[0]=37.0/378.0;
|
|
state.m_rkc[1]=0;
|
|
state.m_rkc[2]=250.0/621.0;
|
|
state.m_rkc[3]=125.0/594.0;
|
|
state.m_rkc[4]=0;
|
|
state.m_rkc[5]=512.0/1771.0;
|
|
//--- allocation
|
|
ArrayResize(state.m_rkcs,6);
|
|
//--- calculation
|
|
state.m_rkcs[0]=2825.0/27648.0;
|
|
state.m_rkcs[1]=0;
|
|
state.m_rkcs[2]=18575.0/48384.0;
|
|
state.m_rkcs[3]=13525.0/55296.0;
|
|
state.m_rkcs[4]=277.0/14336.0;
|
|
state.m_rkcs[5]=1.0/4.0;
|
|
state.m_rkk.Resize(6,n);
|
|
//--- Main cycle consists of two iterations:
|
|
//--- * outer where we travel from X[i-1] to X[i]
|
|
//--- * inner where we travel inside [X[i-1],X[i]]
|
|
state.m_ytbl.Resize(m,n);
|
|
ArrayResize(state.m_escale,n);
|
|
ArrayResize(state.m_yn,n);
|
|
ArrayResize(state.m_yns,n);
|
|
//--- change value
|
|
xc=state.m_xg[0];
|
|
for(i_=0; i_<=n-1; i_++)
|
|
state.m_ytbl.Set(0,i_,state.m_yc[i_]);
|
|
for(j=0; j<=n-1; j++)
|
|
state.m_escale[j]=0;
|
|
i=1;
|
|
//--- check
|
|
if(i>m-1)
|
|
{
|
|
state.m_repterminationtype=1;
|
|
//--- return result
|
|
return(false);
|
|
}
|
|
//--- begin inner iteration
|
|
return(Func_lbl_6(state,n,m,i,j,k,klimit,gridpoint,xc,v,h,h2,err,maxgrowpow));
|
|
}
|
|
//+------------------------------------------------------------------+
|
|
//| Auxiliary function for ODESolverIteration. Is a product to get |
|
|
//| rid of the operator unconditional jump goto |
|
|
//+------------------------------------------------------------------+
|
|
void CODESolver::Func_lbl_rcomm(CODESolverState &state,int n,int m,
|
|
int i,int j,int k,int klimit,
|
|
bool gridpoint,double xc,double v,
|
|
double h,double h2,double err,
|
|
double maxgrowpow)
|
|
{
|
|
//--- save
|
|
state.m_rstate.ia.Set(0,n);
|
|
state.m_rstate.ia.Set(1,m);
|
|
state.m_rstate.ia.Set(2,i);
|
|
state.m_rstate.ia.Set(3,j);
|
|
state.m_rstate.ia.Set(4,k);
|
|
state.m_rstate.ia.Set(5,klimit);
|
|
state.m_rstate.ba[0]= gridpoint;
|
|
state.m_rstate.ra.Set(0,xc);
|
|
state.m_rstate.ra.Set(1,v);
|
|
state.m_rstate.ra.Set(2,h);
|
|
state.m_rstate.ra.Set(3,h2);
|
|
state.m_rstate.ra.Set(4,err);
|
|
state.m_rstate.ra.Set(5,maxgrowpow);
|
|
}
|
|
//+------------------------------------------------------------------+
|
|
//| Auxiliary function for ODESolverIteration. Is a product to get |
|
|
//| rid of the operator unconditional jump goto |
|
|
//+------------------------------------------------------------------+
|
|
bool CODESolver::Func_lbl_6(CODESolverState &state,int &n,int &m,
|
|
int &i,int &j,int &k,int &klimit,
|
|
bool &gridpoint,double &xc,double &v,
|
|
double &h,double &h2,double &err,
|
|
double &maxgrowpow)
|
|
{
|
|
//--- truncate step if needed (beyond right boundary).
|
|
//--- determine should we store X or not
|
|
if(xc+h>=state.m_xg[i])
|
|
{
|
|
h=state.m_xg[i]-xc;
|
|
gridpoint=true;
|
|
}
|
|
else
|
|
gridpoint=false;
|
|
//--- Update error scale maximums
|
|
//--- These maximums are initialized by zeros,
|
|
//--- then updated every iterations.
|
|
for(j=0; j<=n-1; j++)
|
|
state.m_escale[j]=MathMax(state.m_escale[j],MathAbs(state.m_yc[j]));
|
|
//--- make one step:
|
|
//--- 1. calculate all info needed to do step
|
|
//--- 2. update errors scale maximums using values/derivatives
|
|
//--- obtained during (1)
|
|
//--- Take into account that we use scaling of X to reduce task
|
|
//--- to the form where x[0] < x[1] < ... < x[n-1]. So X is
|
|
//--- replaced by x=xscale*t,and dy/dx=f(y,x) is replaced
|
|
//--- by dy/dt=xscale*f(y,xscale*t).
|
|
for(int i_=0; i_<=n-1; i_++)
|
|
state.m_yn[i_]=state.m_yc[i_];
|
|
for(int i_=0; i_<=n-1; i_++)
|
|
state.m_yns[i_]=state.m_yc[i_];
|
|
k=0;
|
|
//--- function call, return result
|
|
return(Func_lbl_8(state,n,m,i,j,k,klimit,gridpoint,xc,v,h,h2,err,maxgrowpow));
|
|
}
|
|
//+------------------------------------------------------------------+
|
|
//| Auxiliary function for ODESolverIteration. Is a product to get |
|
|
//| rid of the operator unconditional jump goto |
|
|
//+------------------------------------------------------------------+
|
|
bool CODESolver::Func_lbl_8(CODESolverState &state,int &n,int &m,
|
|
int &i,int &j,int &k,int &klimit,
|
|
bool &gridpoint,double &xc,double &v,
|
|
double &h,double &h2,double &err,
|
|
double &maxgrowpow)
|
|
{
|
|
//--- check
|
|
if(k>5)
|
|
return(Func_lbl_10(state,n,m,i,j,k,klimit,gridpoint,xc,v,h,h2,err,maxgrowpow));
|
|
//--- prepare data for the next update of YN/YNS
|
|
state.m_x=state.m_xscale*(xc+state.m_rka[k]*h);
|
|
for(int i_=0; i_<=n-1; i_++)
|
|
state.m_y[i_]=state.m_yc[i_];
|
|
//--- calculation
|
|
for(j=0; j<=k-1; j++)
|
|
{
|
|
v=state.m_rkb[k][j];
|
|
for(int i_=0; i_<=n-1; i_++)
|
|
state.m_y[i_]=state.m_y[i_]+v*state.m_rkk[j][i_];
|
|
}
|
|
state.m_needdy=true;
|
|
state.m_rstate.stage=0;
|
|
//--- Saving state
|
|
Func_lbl_rcomm(state,n,m,i,j,k,klimit,gridpoint,xc,v,h,h2,err,maxgrowpow);
|
|
//--- return result
|
|
return(true);
|
|
}
|
|
//+------------------------------------------------------------------+
|
|
//| Auxiliary function for ODESolverIteration. Is a product to get |
|
|
//| rid of the operator unconditional jump goto |
|
|
//+------------------------------------------------------------------+
|
|
bool CODESolver::Func_lbl_10(CODESolverState &state,int &n,int &m,
|
|
int &i,int &j,int &k,int &klimit,
|
|
bool &gridpoint,double &xc,double &v,
|
|
double &h,double &h2,double &err,
|
|
double &maxgrowpow)
|
|
{
|
|
//--- estimate error
|
|
err=0;
|
|
for(j=0; j<=n-1; j++)
|
|
{
|
|
//--- check
|
|
if(!state.m_fraceps)
|
|
{
|
|
//--- absolute error is estimated
|
|
err=MathMax(err,MathAbs(state.m_yn[j]-state.m_yns[j]));
|
|
}
|
|
else
|
|
{
|
|
//--- Relative error is estimated
|
|
v=state.m_escale[j];
|
|
//--- check
|
|
if(v==0.0)
|
|
v=1;
|
|
err=MathMax(err,MathAbs(state.m_yn[j]-state.m_yns[j])/v);
|
|
}
|
|
}
|
|
//--- calculate new step,restart if necessary
|
|
if(maxgrowpow*err<=state.m_eps)
|
|
h2=m_odesolvermaxgrow*h;
|
|
else
|
|
h2=h*MathPow(state.m_eps/err,0.2);
|
|
//--- check
|
|
if(h2<h/m_odesolvermaxshrink)
|
|
h2=h/m_odesolvermaxshrink;
|
|
//--- check
|
|
if(err>state.m_eps)
|
|
{
|
|
h=h2;
|
|
//--- begin inner iteration
|
|
return(Func_lbl_6(state,n,m,i,j,k,klimit,gridpoint,xc,v,h,h2,err,maxgrowpow));
|
|
}
|
|
//--- advance position
|
|
xc=xc+h;
|
|
for(int i_=0; i_<=n-1; i_++)
|
|
state.m_yc[i_]=state.m_yn[i_];
|
|
//--- update H
|
|
h=h2;
|
|
//--- break on grid point
|
|
if(gridpoint)
|
|
{
|
|
//--- save result
|
|
for(int i_=0; i_<=n-1; i_++)
|
|
state.m_ytbl.Set(i,i_,state.m_yc[i_]);
|
|
i=i+1;
|
|
//--- check
|
|
if(i>m-1)
|
|
{
|
|
state.m_repterminationtype=1;
|
|
return(false);
|
|
}
|
|
//--- begin inner iteration
|
|
return(Func_lbl_6(state,n,m,i,j,k,klimit,gridpoint,xc,v,h,h2,err,maxgrowpow));
|
|
}
|
|
//--- begin inner iteration
|
|
return(Func_lbl_6(state,n,m,i,j,k,klimit,gridpoint,xc,v,h,h2,err,maxgrowpow));
|
|
}
|
|
//+------------------------------------------------------------------+
|