1 /************************************************************************* 2 ALGLIB 3.15.0 (source code generated 2019-02-20) 3 Copyright (c) Sergey Bochkanov (ALGLIB project). 4 5 >>> SOURCE LICENSE >>> 6 This program is free software; you can redistribute it and/or modify 7 it under the terms of the GNU General Public License as published by 8 the Free Software Foundation (www.fsf.org); either version 2 of the 9 License, or (at your option) any later version. 10 11 This program is distributed in the hope that it will be useful, 12 but WITHOUT ANY WARRANTY; without even the implied warranty of 13 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 14 GNU General Public License for more details. 15 16 A copy of the GNU General Public License is available at 17 http://www.fsf.org/licensing/licenses 18 >>> END OF LICENSE >>> 19 *************************************************************************/ 20 #ifndef _diffequations_pkg_h 21 #define _diffequations_pkg_h 22 #include "ap.h" 23 #include "alglibinternal.h" 24 25 ///////////////////////////////////////////////////////////////////////// 26 // 27 // THIS SECTION CONTAINS COMPUTATIONAL CORE DECLARATIONS (DATATYPES) 28 // 29 ///////////////////////////////////////////////////////////////////////// 30 namespace alglib_impl 31 { 32 #if defined(AE_COMPILE_ODESOLVER) || !defined(AE_PARTIAL_BUILD) 33 typedef struct 34 { 35 ae_int_t n; 36 ae_int_t m; 37 double xscale; 38 double h; 39 double eps; 40 ae_bool fraceps; 41 ae_vector yc; 42 ae_vector escale; 43 ae_vector xg; 44 ae_int_t solvertype; 45 ae_bool needdy; 46 double x; 47 ae_vector y; 48 ae_vector dy; 49 ae_matrix ytbl; 50 ae_int_t repterminationtype; 51 ae_int_t repnfev; 52 ae_vector yn; 53 ae_vector yns; 54 ae_vector rka; 55 ae_vector rkc; 56 ae_vector rkcs; 57 ae_matrix rkb; 58 ae_matrix rkk; 59 rcommstate rstate; 60 } odesolverstate; 61 typedef struct 62 { 63 ae_int_t nfev; 64 ae_int_t terminationtype; 65 } odesolverreport; 66 #endif 67 68 } 69 70 ///////////////////////////////////////////////////////////////////////// 71 // 72 // THIS SECTION CONTAINS C++ INTERFACE 73 // 74 ///////////////////////////////////////////////////////////////////////// 75 namespace alglib 76 { 77 78 #if defined(AE_COMPILE_ODESOLVER) || !defined(AE_PARTIAL_BUILD) 79 /************************************************************************* 80 81 *************************************************************************/ 82 class _odesolverstate_owner 83 { 84 public: 85 _odesolverstate_owner(); 86 _odesolverstate_owner(const _odesolverstate_owner &rhs); 87 _odesolverstate_owner& operator=(const _odesolverstate_owner &rhs); 88 virtual ~_odesolverstate_owner(); 89 alglib_impl::odesolverstate* c_ptr(); 90 alglib_impl::odesolverstate* c_ptr() const; 91 protected: 92 alglib_impl::odesolverstate *p_struct; 93 }; 94 class odesolverstate : public _odesolverstate_owner 95 { 96 public: 97 odesolverstate(); 98 odesolverstate(const odesolverstate &rhs); 99 odesolverstate& operator=(const odesolverstate &rhs); 100 virtual ~odesolverstate(); 101 ae_bool &needdy; 102 real_1d_array y; 103 real_1d_array dy; 104 double &x; 105 106 }; 107 108 109 /************************************************************************* 110 111 *************************************************************************/ 112 class _odesolverreport_owner 113 { 114 public: 115 _odesolverreport_owner(); 116 _odesolverreport_owner(const _odesolverreport_owner &rhs); 117 _odesolverreport_owner& operator=(const _odesolverreport_owner &rhs); 118 virtual ~_odesolverreport_owner(); 119 alglib_impl::odesolverreport* c_ptr(); 120 alglib_impl::odesolverreport* c_ptr() const; 121 protected: 122 alglib_impl::odesolverreport *p_struct; 123 }; 124 class odesolverreport : public _odesolverreport_owner 125 { 126 public: 127 odesolverreport(); 128 odesolverreport(const odesolverreport &rhs); 129 odesolverreport& operator=(const odesolverreport &rhs); 130 virtual ~odesolverreport(); 131 ae_int_t &nfev; 132 ae_int_t &terminationtype; 133 134 }; 135 #endif 136 137 #if defined(AE_COMPILE_ODESOLVER) || !defined(AE_PARTIAL_BUILD) 138 /************************************************************************* 139 Cash-Karp adaptive ODE solver. 140 141 This subroutine solves ODE Y'=f(Y,x) with initial conditions Y(xs)=Ys 142 (here Y may be single variable or vector of N variables). 143 144 INPUT PARAMETERS: 145 Y - initial conditions, array[0..N-1]. 146 contains values of Y[] at X[0] 147 N - system size 148 X - points at which Y should be tabulated, array[0..M-1] 149 integrations starts at X[0], ends at X[M-1], intermediate 150 values at X[i] are returned too. 151 SHOULD BE ORDERED BY ASCENDING OR BY DESCENDING! 152 M - number of intermediate points + first point + last point: 153 * M>2 means that you need both Y(X[M-1]) and M-2 values at 154 intermediate points 155 * M=2 means that you want just to integrate from X[0] to 156 X[1] and don't interested in intermediate values. 157 * M=1 means that you don't want to integrate :) 158 it is degenerate case, but it will be handled correctly. 159 * M<1 means error 160 Eps - tolerance (absolute/relative error on each step will be 161 less than Eps). When passing: 162 * Eps>0, it means desired ABSOLUTE error 163 * Eps<0, it means desired RELATIVE error. Relative errors 164 are calculated with respect to maximum values of Y seen 165 so far. Be careful to use this criterion when starting 166 from Y[] that are close to zero. 167 H - initial step lenth, it will be adjusted automatically 168 after the first step. If H=0, step will be selected 169 automatically (usualy it will be equal to 0.001 of 170 min(x[i]-x[j])). 171 172 OUTPUT PARAMETERS 173 State - structure which stores algorithm state between subsequent 174 calls of OdeSolverIteration. Used for reverse communication. 175 This structure should be passed to the OdeSolverIteration 176 subroutine. 177 178 SEE ALSO 179 AutoGKSmoothW, AutoGKSingular, AutoGKIteration, AutoGKResults. 180 181 182 -- ALGLIB -- 183 Copyright 01.09.2009 by Bochkanov Sergey 184 *************************************************************************/ 185 void odesolverrkck(const real_1d_array &y, const ae_int_t n, const real_1d_array &x, const ae_int_t m, const double eps, const double h, odesolverstate &state, const xparams _xparams = alglib::xdefault); 186 void odesolverrkck(const real_1d_array &y, const real_1d_array &x, const double eps, const double h, odesolverstate &state, const xparams _xparams = alglib::xdefault); 187 188 189 /************************************************************************* 190 This function provides reverse communication interface 191 Reverse communication interface is not documented or recommended to use. 192 See below for functions which provide better documented API 193 *************************************************************************/ 194 bool odesolveriteration(const odesolverstate &state, const xparams _xparams = alglib::xdefault); 195 196 197 /************************************************************************* 198 This function is used to launcn iterations of ODE solver 199 200 It accepts following parameters: 201 diff - callback which calculates dy/dx for given y and x 202 ptr - optional pointer which is passed to diff; can be NULL 203 204 205 -- ALGLIB -- 206 Copyright 01.09.2009 by Bochkanov Sergey 207 208 *************************************************************************/ 209 void odesolversolve(odesolverstate &state, 210 void (*diff)(const real_1d_array &y, double x, real_1d_array &dy, void *ptr), 211 void *ptr = NULL, const xparams _xparams = alglib::xdefault); 212 213 214 /************************************************************************* 215 ODE solver results 216 217 Called after OdeSolverIteration returned False. 218 219 INPUT PARAMETERS: 220 State - algorithm state (used by OdeSolverIteration). 221 222 OUTPUT PARAMETERS: 223 M - number of tabulated values, M>=1 224 XTbl - array[0..M-1], values of X 225 YTbl - array[0..M-1,0..N-1], values of Y in X[i] 226 Rep - solver report: 227 * Rep.TerminationType completetion code: 228 * -2 X is not ordered by ascending/descending or 229 there are non-distinct X[], i.e. X[i]=X[i+1] 230 * -1 incorrect parameters were specified 231 * 1 task has been solved 232 * Rep.NFEV contains number of function calculations 233 234 -- ALGLIB -- 235 Copyright 01.09.2009 by Bochkanov Sergey 236 *************************************************************************/ 237 void odesolverresults(const odesolverstate &state, ae_int_t &m, real_1d_array &xtbl, real_2d_array &ytbl, odesolverreport &rep, const xparams _xparams = alglib::xdefault); 238 #endif 239 } 240 241 ///////////////////////////////////////////////////////////////////////// 242 // 243 // THIS SECTION CONTAINS COMPUTATIONAL CORE DECLARATIONS (FUNCTIONS) 244 // 245 ///////////////////////////////////////////////////////////////////////// 246 namespace alglib_impl 247 { 248 #if defined(AE_COMPILE_ODESOLVER) || !defined(AE_PARTIAL_BUILD) 249 void odesolverrkck(/* Real */ ae_vector* y, 250 ae_int_t n, 251 /* Real */ ae_vector* x, 252 ae_int_t m, 253 double eps, 254 double h, 255 odesolverstate* state, 256 ae_state *_state); 257 ae_bool odesolveriteration(odesolverstate* state, ae_state *_state); 258 void odesolverresults(odesolverstate* state, 259 ae_int_t* m, 260 /* Real */ ae_vector* xtbl, 261 /* Real */ ae_matrix* ytbl, 262 odesolverreport* rep, 263 ae_state *_state); 264 void _odesolverstate_init(void* _p, ae_state *_state, ae_bool make_automatic); 265 void _odesolverstate_init_copy(void* _dst, void* _src, ae_state *_state, ae_bool make_automatic); 266 void _odesolverstate_clear(void* _p); 267 void _odesolverstate_destroy(void* _p); 268 void _odesolverreport_init(void* _p, ae_state *_state, ae_bool make_automatic); 269 void _odesolverreport_init_copy(void* _dst, void* _src, ae_state *_state, ae_bool make_automatic); 270 void _odesolverreport_clear(void* _p); 271 void _odesolverreport_destroy(void* _p); 272 #endif 273 274 } 275 #endif 276 277