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