Princeec13/mql5
Watch
0
0
Fork
You've already forked mql5
 1
Code Issues Pull requests 1 Projects Releases Packages Wiki Activity Actions
mql5/Include/Math/Alglib/diffequations.mqh

849 lines
32 KiB
MQL5
RawPermalinkNormal ViewHistory

new files added
2025-08-05 01:57:33 -04:00
//+------------------------------------------------------------------+
//| 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));
}
//+------------------------------------------------------------------+
Reference in a new issueCopy permalink