//+------------------------------------------------------------------+ //| 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)=m,"ODESolverRKCK: Length(X)=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[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(h2state.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)); } //+------------------------------------------------------------------+