Free cookie consent management tool by TermsFeed Policy Generator

source: trunk/sources/HeuristicLab.ExtLibs/HeuristicLab.ALGLIB/3.5.0/ALGLIB-3.5.0/diffequations.cs @ 7673

Last change on this file since 7673 was 7673, checked in by mkommend, 11 years ago

#1808: Added new alglib release 3.5.0 to the external libraries.

File size: 32.7 KB
Line 
1/*************************************************************************
2Copyright (c) Sergey Bochkanov (ALGLIB project).
3
4>>> SOURCE LICENSE >>>
5This program is free software; you can redistribute it and/or modify
6it under the terms of the GNU General Public License as published by
7the Free Software Foundation (www.fsf.org); either version 2 of the
8License, or (at your option) any later version.
9
10This program is distributed in the hope that it will be useful,
11but WITHOUT ANY WARRANTY; without even the implied warranty of
12MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
13GNU General Public License for more details.
14
15A copy of the GNU General Public License is available at
16http://www.fsf.org/licensing/licenses
17>>> END OF LICENSE >>>
18*************************************************************************/
19#pragma warning disable 162
20#pragma warning disable 219
21using System;
22
23public partial class alglib
24{
25
26
27    /*************************************************************************
28
29    *************************************************************************/
30    public class odesolverstate
31    {
32        //
33        // Public declarations
34        //
35        public bool needdy { get { return _innerobj.needdy; } set { _innerobj.needdy = value; } }
36        public double[] y { get { return _innerobj.y; } }
37        public double[] dy { get { return _innerobj.dy; } }
38        public double x { get { return _innerobj.x; } set { _innerobj.x = value; } }
39
40        public odesolverstate()
41        {
42            _innerobj = new odesolver.odesolverstate();
43        }
44
45        //
46        // Although some of declarations below are public, you should not use them
47        // They are intended for internal use only
48        //
49        private odesolver.odesolverstate _innerobj;
50        public odesolver.odesolverstate innerobj { get { return _innerobj; } }
51        public odesolverstate(odesolver.odesolverstate obj)
52        {
53            _innerobj = obj;
54        }
55    }
56
57
58    /*************************************************************************
59
60    *************************************************************************/
61    public class odesolverreport
62    {
63        //
64        // Public declarations
65        //
66        public int nfev { get { return _innerobj.nfev; } set { _innerobj.nfev = value; } }
67        public int terminationtype { get { return _innerobj.terminationtype; } set { _innerobj.terminationtype = value; } }
68
69        public odesolverreport()
70        {
71            _innerobj = new odesolver.odesolverreport();
72        }
73
74        //
75        // Although some of declarations below are public, you should not use them
76        // They are intended for internal use only
77        //
78        private odesolver.odesolverreport _innerobj;
79        public odesolver.odesolverreport innerobj { get { return _innerobj; } }
80        public odesolverreport(odesolver.odesolverreport obj)
81        {
82            _innerobj = obj;
83        }
84    }
85
86    /*************************************************************************
87    Cash-Karp adaptive ODE solver.
88
89    This subroutine solves ODE  Y'=f(Y,x)  with  initial  conditions  Y(xs)=Ys
90    (here Y may be single variable or vector of N variables).
91
92    INPUT PARAMETERS:
93        Y       -   initial conditions, array[0..N-1].
94                    contains values of Y[] at X[0]
95        N       -   system size
96        X       -   points at which Y should be tabulated, array[0..M-1]
97                    integrations starts at X[0], ends at X[M-1],  intermediate
98                    values at X[i] are returned too.
99                    SHOULD BE ORDERED BY ASCENDING OR BY DESCENDING!!!!
100        M       -   number of intermediate points + first point + last point:
101                    * M>2 means that you need both Y(X[M-1]) and M-2 values at
102                      intermediate points
103                    * M=2 means that you want just to integrate from  X[0]  to
104                      X[1] and don't interested in intermediate values.
105                    * M=1 means that you don't want to integrate :)
106                      it is degenerate case, but it will be handled correctly.
107                    * M<1 means error
108        Eps     -   tolerance (absolute/relative error on each  step  will  be
109                    less than Eps). When passing:
110                    * Eps>0, it means desired ABSOLUTE error
111                    * Eps<0, it means desired RELATIVE error.  Relative errors
112                      are calculated with respect to maximum values of  Y seen
113                      so far. Be careful to use this criterion  when  starting
114                      from Y[] that are close to zero.
115        H       -   initial  step  lenth,  it  will  be adjusted automatically
116                    after the first  step.  If  H=0,  step  will  be  selected
117                    automatically  (usualy  it  will  be  equal  to  0.001  of
118                    min(x[i]-x[j])).
119
120    OUTPUT PARAMETERS
121        State   -   structure which stores algorithm state between  subsequent
122                    calls of OdeSolverIteration. Used for reverse communication.
123                    This structure should be passed  to the OdeSolverIteration
124                    subroutine.
125
126    SEE ALSO
127        AutoGKSmoothW, AutoGKSingular, AutoGKIteration, AutoGKResults.
128
129
130      -- ALGLIB --
131         Copyright 01.09.2009 by Bochkanov Sergey
132    *************************************************************************/
133    public static void odesolverrkck(double[] y, int n, double[] x, int m, double eps, double h, out odesolverstate state)
134    {
135        state = new odesolverstate();
136        odesolver.odesolverrkck(y, n, x, m, eps, h, state.innerobj);
137        return;
138    }
139    public static void odesolverrkck(double[] y, double[] x, double eps, double h, out odesolverstate state)
140    {
141        int n;
142        int m;
143
144        state = new odesolverstate();
145        n = ap.len(y);
146        m = ap.len(x);
147        odesolver.odesolverrkck(y, n, x, m, eps, h, state.innerobj);
148
149        return;
150    }
151
152    /*************************************************************************
153    This function provides reverse communication interface
154    Reverse communication interface is not documented or recommended to use.
155    See below for functions which provide better documented API
156    *************************************************************************/
157    public static bool odesolveriteration(odesolverstate state)
158    {
159
160        bool result = odesolver.odesolveriteration(state.innerobj);
161        return result;
162    }
163    /*************************************************************************
164    This function is used to launcn iterations of ODE solver
165
166    It accepts following parameters:
167        diff    -   callback which calculates dy/dx for given y and x
168        obj     -   optional object which is passed to diff; can be NULL
169
170
171      -- ALGLIB --
172         Copyright 01.09.2009 by Bochkanov Sergey
173
174    *************************************************************************/
175    public static void odesolversolve(odesolverstate state, ndimensional_ode_rp diff, object obj)
176    {
177        if( diff==null )
178            throw new alglibexception("ALGLIB: error in 'odesolversolve()' (diff is null)");
179        while( alglib.odesolveriteration(state) )
180        {
181            if( state.needdy )
182            {
183                diff(state.innerobj.y, state.innerobj.x, state.innerobj.dy, obj);
184                continue;
185            }
186            throw new alglibexception("ALGLIB: unexpected error in 'odesolversolve'");
187        }
188    }
189
190
191
192    /*************************************************************************
193    ODE solver results
194
195    Called after OdeSolverIteration returned False.
196
197    INPUT PARAMETERS:
198        State   -   algorithm state (used by OdeSolverIteration).
199
200    OUTPUT PARAMETERS:
201        M       -   number of tabulated values, M>=1
202        XTbl    -   array[0..M-1], values of X
203        YTbl    -   array[0..M-1,0..N-1], values of Y in X[i]
204        Rep     -   solver report:
205                    * Rep.TerminationType completetion code:
206                        * -2    X is not ordered  by  ascending/descending  or
207                                there are non-distinct X[],  i.e.  X[i]=X[i+1]
208                        * -1    incorrect parameters were specified
209                        *  1    task has been solved
210                    * Rep.NFEV contains number of function calculations
211
212      -- ALGLIB --
213         Copyright 01.09.2009 by Bochkanov Sergey
214    *************************************************************************/
215    public static void odesolverresults(odesolverstate state, out int m, out double[] xtbl, out double[,] ytbl, out odesolverreport rep)
216    {
217        m = 0;
218        xtbl = new double[0];
219        ytbl = new double[0,0];
220        rep = new odesolverreport();
221        odesolver.odesolverresults(state.innerobj, ref m, ref xtbl, ref ytbl, rep.innerobj);
222        return;
223    }
224
225}
226public partial class alglib
227{
228    public class odesolver
229    {
230        public class odesolverstate
231        {
232            public int n;
233            public int m;
234            public double xscale;
235            public double h;
236            public double eps;
237            public bool fraceps;
238            public double[] yc;
239            public double[] escale;
240            public double[] xg;
241            public int solvertype;
242            public bool needdy;
243            public double x;
244            public double[] y;
245            public double[] dy;
246            public double[,] ytbl;
247            public int repterminationtype;
248            public int repnfev;
249            public double[] yn;
250            public double[] yns;
251            public double[] rka;
252            public double[] rkc;
253            public double[] rkcs;
254            public double[,] rkb;
255            public double[,] rkk;
256            public rcommstate rstate;
257            public odesolverstate()
258            {
259                yc = new double[0];
260                escale = new double[0];
261                xg = new double[0];
262                y = new double[0];
263                dy = new double[0];
264                ytbl = new double[0,0];
265                yn = new double[0];
266                yns = new double[0];
267                rka = new double[0];
268                rkc = new double[0];
269                rkcs = new double[0];
270                rkb = new double[0,0];
271                rkk = new double[0,0];
272                rstate = new rcommstate();
273            }
274        };
275
276
277        public class odesolverreport
278        {
279            public int nfev;
280            public int terminationtype;
281        };
282
283
284
285
286        public const double odesolvermaxgrow = 3.0;
287        public const double odesolvermaxshrink = 10.0;
288
289
290        /*************************************************************************
291        Cash-Karp adaptive ODE solver.
292
293        This subroutine solves ODE  Y'=f(Y,x)  with  initial  conditions  Y(xs)=Ys
294        (here Y may be single variable or vector of N variables).
295
296        INPUT PARAMETERS:
297            Y       -   initial conditions, array[0..N-1].
298                        contains values of Y[] at X[0]
299            N       -   system size
300            X       -   points at which Y should be tabulated, array[0..M-1]
301                        integrations starts at X[0], ends at X[M-1],  intermediate
302                        values at X[i] are returned too.
303                        SHOULD BE ORDERED BY ASCENDING OR BY DESCENDING!!!!
304            M       -   number of intermediate points + first point + last point:
305                        * M>2 means that you need both Y(X[M-1]) and M-2 values at
306                          intermediate points
307                        * M=2 means that you want just to integrate from  X[0]  to
308                          X[1] and don't interested in intermediate values.
309                        * M=1 means that you don't want to integrate :)
310                          it is degenerate case, but it will be handled correctly.
311                        * M<1 means error
312            Eps     -   tolerance (absolute/relative error on each  step  will  be
313                        less than Eps). When passing:
314                        * Eps>0, it means desired ABSOLUTE error
315                        * Eps<0, it means desired RELATIVE error.  Relative errors
316                          are calculated with respect to maximum values of  Y seen
317                          so far. Be careful to use this criterion  when  starting
318                          from Y[] that are close to zero.
319            H       -   initial  step  lenth,  it  will  be adjusted automatically
320                        after the first  step.  If  H=0,  step  will  be  selected
321                        automatically  (usualy  it  will  be  equal  to  0.001  of
322                        min(x[i]-x[j])).
323
324        OUTPUT PARAMETERS
325            State   -   structure which stores algorithm state between  subsequent
326                        calls of OdeSolverIteration. Used for reverse communication.
327                        This structure should be passed  to the OdeSolverIteration
328                        subroutine.
329
330        SEE ALSO
331            AutoGKSmoothW, AutoGKSingular, AutoGKIteration, AutoGKResults.
332
333
334          -- ALGLIB --
335             Copyright 01.09.2009 by Bochkanov Sergey
336        *************************************************************************/
337        public static void odesolverrkck(double[] y,
338            int n,
339            double[] x,
340            int m,
341            double eps,
342            double h,
343            odesolverstate state)
344        {
345            alglib.ap.assert(n>=1, "ODESolverRKCK: N<1!");
346            alglib.ap.assert(m>=1, "ODESolverRKCK: M<1!");
347            alglib.ap.assert(alglib.ap.len(y)>=n, "ODESolverRKCK: Length(Y)<N!");
348            alglib.ap.assert(alglib.ap.len(x)>=m, "ODESolverRKCK: Length(X)<M!");
349            alglib.ap.assert(apserv.isfinitevector(y, n), "ODESolverRKCK: Y contains infinite or NaN values!");
350            alglib.ap.assert(apserv.isfinitevector(x, m), "ODESolverRKCK: Y contains infinite or NaN values!");
351            alglib.ap.assert(math.isfinite(eps), "ODESolverRKCK: Eps is not finite!");
352            alglib.ap.assert((double)(eps)!=(double)(0), "ODESolverRKCK: Eps is zero!");
353            alglib.ap.assert(math.isfinite(h), "ODESolverRKCK: H is not finite!");
354            odesolverinit(0, y, n, x, m, eps, h, state);
355        }
356
357
358        /*************************************************************************
359
360          -- ALGLIB --
361             Copyright 01.09.2009 by Bochkanov Sergey
362        *************************************************************************/
363        public static bool odesolveriteration(odesolverstate state)
364        {
365            bool result = new bool();
366            int n = 0;
367            int m = 0;
368            int i = 0;
369            int j = 0;
370            int k = 0;
371            double xc = 0;
372            double v = 0;
373            double h = 0;
374            double h2 = 0;
375            bool gridpoint = new bool();
376            double err = 0;
377            double maxgrowpow = 0;
378            int klimit = 0;
379            int i_ = 0;
380
381           
382            //
383            // Reverse communication preparations
384            // I know it looks ugly, but it works the same way
385            // anywhere from C++ to Python.
386            //
387            // This code initializes locals by:
388            // * random values determined during code
389            //   generation - on first subroutine call
390            // * values from previous call - on subsequent calls
391            //
392            if( state.rstate.stage>=0 )
393            {
394                n = state.rstate.ia[0];
395                m = state.rstate.ia[1];
396                i = state.rstate.ia[2];
397                j = state.rstate.ia[3];
398                k = state.rstate.ia[4];
399                klimit = state.rstate.ia[5];
400                gridpoint = state.rstate.ba[0];
401                xc = state.rstate.ra[0];
402                v = state.rstate.ra[1];
403                h = state.rstate.ra[2];
404                h2 = state.rstate.ra[3];
405                err = state.rstate.ra[4];
406                maxgrowpow = state.rstate.ra[5];
407            }
408            else
409            {
410                n = -983;
411                m = -989;
412                i = -834;
413                j = 900;
414                k = -287;
415                klimit = 364;
416                gridpoint = false;
417                xc = -338;
418                v = -686;
419                h = 912;
420                h2 = 585;
421                err = 497;
422                maxgrowpow = -271;
423            }
424            if( state.rstate.stage==0 )
425            {
426                goto lbl_0;
427            }
428           
429            //
430            // Routine body
431            //
432           
433            //
434            // prepare
435            //
436            if( state.repterminationtype!=0 )
437            {
438                result = false;
439                return result;
440            }
441            n = state.n;
442            m = state.m;
443            h = state.h;
444            maxgrowpow = Math.Pow(odesolvermaxgrow, 5);
445            state.repnfev = 0;
446           
447            //
448            // some preliminary checks for internal errors
449            // after this we assume that H>0 and M>1
450            //
451            alglib.ap.assert((double)(state.h)>(double)(0), "ODESolver: internal error");
452            alglib.ap.assert(m>1, "ODESolverIteration: internal error");
453           
454            //
455            // choose solver
456            //
457            if( state.solvertype!=0 )
458            {
459                goto lbl_1;
460            }
461           
462            //
463            // Cask-Karp solver
464            // Prepare coefficients table.
465            // Check it for errors
466            //
467            state.rka = new double[6];
468            state.rka[0] = 0;
469            state.rka[1] = (double)1/(double)5;
470            state.rka[2] = (double)3/(double)10;
471            state.rka[3] = (double)3/(double)5;
472            state.rka[4] = 1;
473            state.rka[5] = (double)7/(double)8;
474            state.rkb = new double[6, 5];
475            state.rkb[1,0] = (double)1/(double)5;
476            state.rkb[2,0] = (double)3/(double)40;
477            state.rkb[2,1] = (double)9/(double)40;
478            state.rkb[3,0] = (double)3/(double)10;
479            state.rkb[3,1] = -((double)9/(double)10);
480            state.rkb[3,2] = (double)6/(double)5;
481            state.rkb[4,0] = -((double)11/(double)54);
482            state.rkb[4,1] = (double)5/(double)2;
483            state.rkb[4,2] = -((double)70/(double)27);
484            state.rkb[4,3] = (double)35/(double)27;
485            state.rkb[5,0] = (double)1631/(double)55296;
486            state.rkb[5,1] = (double)175/(double)512;
487            state.rkb[5,2] = (double)575/(double)13824;
488            state.rkb[5,3] = (double)44275/(double)110592;
489            state.rkb[5,4] = (double)253/(double)4096;
490            state.rkc = new double[6];
491            state.rkc[0] = (double)37/(double)378;
492            state.rkc[1] = 0;
493            state.rkc[2] = (double)250/(double)621;
494            state.rkc[3] = (double)125/(double)594;
495            state.rkc[4] = 0;
496            state.rkc[5] = (double)512/(double)1771;
497            state.rkcs = new double[6];
498            state.rkcs[0] = (double)2825/(double)27648;
499            state.rkcs[1] = 0;
500            state.rkcs[2] = (double)18575/(double)48384;
501            state.rkcs[3] = (double)13525/(double)55296;
502            state.rkcs[4] = (double)277/(double)14336;
503            state.rkcs[5] = (double)1/(double)4;
504            state.rkk = new double[6, n];
505           
506            //
507            // Main cycle consists of two iterations:
508            // * outer where we travel from X[i-1] to X[i]
509            // * inner where we travel inside [X[i-1],X[i]]
510            //
511            state.ytbl = new double[m, n];
512            state.escale = new double[n];
513            state.yn = new double[n];
514            state.yns = new double[n];
515            xc = state.xg[0];
516            for(i_=0; i_<=n-1;i_++)
517            {
518                state.ytbl[0,i_] = state.yc[i_];
519            }
520            for(j=0; j<=n-1; j++)
521            {
522                state.escale[j] = 0;
523            }
524            i = 1;
525        lbl_3:
526            if( i>m-1 )
527            {
528                goto lbl_5;
529            }
530           
531            //
532            // begin inner iteration
533            //
534        lbl_6:
535            if( false )
536            {
537                goto lbl_7;
538            }
539           
540            //
541            // truncate step if needed (beyond right boundary).
542            // determine should we store X or not
543            //
544            if( (double)(xc+h)>=(double)(state.xg[i]) )
545            {
546                h = state.xg[i]-xc;
547                gridpoint = true;
548            }
549            else
550            {
551                gridpoint = false;
552            }
553           
554            //
555            // Update error scale maximums
556            //
557            // These maximums are initialized by zeros,
558            // then updated every iterations.
559            //
560            for(j=0; j<=n-1; j++)
561            {
562                state.escale[j] = Math.Max(state.escale[j], Math.Abs(state.yc[j]));
563            }
564           
565            //
566            // make one step:
567            // 1. calculate all info needed to do step
568            // 2. update errors scale maximums using values/derivatives
569            //    obtained during (1)
570            //
571            // Take into account that we use scaling of X to reduce task
572            // to the form where x[0] < x[1] < ... < x[n-1]. So X is
573            // replaced by x=xscale*t, and dy/dx=f(y,x) is replaced
574            // by dy/dt=xscale*f(y,xscale*t).
575            //
576            for(i_=0; i_<=n-1;i_++)
577            {
578                state.yn[i_] = state.yc[i_];
579            }
580            for(i_=0; i_<=n-1;i_++)
581            {
582                state.yns[i_] = state.yc[i_];
583            }
584            k = 0;
585        lbl_8:
586            if( k>5 )
587            {
588                goto lbl_10;
589            }
590           
591            //
592            // prepare data for the next update of YN/YNS
593            //
594            state.x = state.xscale*(xc+state.rka[k]*h);
595            for(i_=0; i_<=n-1;i_++)
596            {
597                state.y[i_] = state.yc[i_];
598            }
599            for(j=0; j<=k-1; j++)
600            {
601                v = state.rkb[k,j];
602                for(i_=0; i_<=n-1;i_++)
603                {
604                    state.y[i_] = state.y[i_] + v*state.rkk[j,i_];
605                }
606            }
607            state.needdy = true;
608            state.rstate.stage = 0;
609            goto lbl_rcomm;
610        lbl_0:
611            state.needdy = false;
612            state.repnfev = state.repnfev+1;
613            v = h*state.xscale;
614            for(i_=0; i_<=n-1;i_++)
615            {
616                state.rkk[k,i_] = v*state.dy[i_];
617            }
618           
619            //
620            // update YN/YNS
621            //
622            v = state.rkc[k];
623            for(i_=0; i_<=n-1;i_++)
624            {
625                state.yn[i_] = state.yn[i_] + v*state.rkk[k,i_];
626            }
627            v = state.rkcs[k];
628            for(i_=0; i_<=n-1;i_++)
629            {
630                state.yns[i_] = state.yns[i_] + v*state.rkk[k,i_];
631            }
632            k = k+1;
633            goto lbl_8;
634        lbl_10:
635           
636            //
637            // estimate error
638            //
639            err = 0;
640            for(j=0; j<=n-1; j++)
641            {
642                if( !state.fraceps )
643                {
644                   
645                    //
646                    // absolute error is estimated
647                    //
648                    err = Math.Max(err, Math.Abs(state.yn[j]-state.yns[j]));
649                }
650                else
651                {
652                   
653                    //
654                    // Relative error is estimated
655                    //
656                    v = state.escale[j];
657                    if( (double)(v)==(double)(0) )
658                    {
659                        v = 1;
660                    }
661                    err = Math.Max(err, Math.Abs(state.yn[j]-state.yns[j])/v);
662                }
663            }
664           
665            //
666            // calculate new step, restart if necessary
667            //
668            if( (double)(maxgrowpow*err)<=(double)(state.eps) )
669            {
670                h2 = odesolvermaxgrow*h;
671            }
672            else
673            {
674                h2 = h*Math.Pow(state.eps/err, 0.2);
675            }
676            if( (double)(h2)<(double)(h/odesolvermaxshrink) )
677            {
678                h2 = h/odesolvermaxshrink;
679            }
680            if( (double)(err)>(double)(state.eps) )
681            {
682                h = h2;
683                goto lbl_6;
684            }
685           
686            //
687            // advance position
688            //
689            xc = xc+h;
690            for(i_=0; i_<=n-1;i_++)
691            {
692                state.yc[i_] = state.yn[i_];
693            }
694           
695            //
696            // update H
697            //
698            h = h2;
699           
700            //
701            // break on grid point
702            //
703            if( gridpoint )
704            {
705                goto lbl_7;
706            }
707            goto lbl_6;
708        lbl_7:
709           
710            //
711            // save result
712            //
713            for(i_=0; i_<=n-1;i_++)
714            {
715                state.ytbl[i,i_] = state.yc[i_];
716            }
717            i = i+1;
718            goto lbl_3;
719        lbl_5:
720            state.repterminationtype = 1;
721            result = false;
722            return result;
723        lbl_1:
724            result = false;
725            return result;
726           
727            //
728            // Saving state
729            //
730        lbl_rcomm:
731            result = true;
732            state.rstate.ia[0] = n;
733            state.rstate.ia[1] = m;
734            state.rstate.ia[2] = i;
735            state.rstate.ia[3] = j;
736            state.rstate.ia[4] = k;
737            state.rstate.ia[5] = klimit;
738            state.rstate.ba[0] = gridpoint;
739            state.rstate.ra[0] = xc;
740            state.rstate.ra[1] = v;
741            state.rstate.ra[2] = h;
742            state.rstate.ra[3] = h2;
743            state.rstate.ra[4] = err;
744            state.rstate.ra[5] = maxgrowpow;
745            return result;
746        }
747
748
749        /*************************************************************************
750        ODE solver results
751
752        Called after OdeSolverIteration returned False.
753
754        INPUT PARAMETERS:
755            State   -   algorithm state (used by OdeSolverIteration).
756
757        OUTPUT PARAMETERS:
758            M       -   number of tabulated values, M>=1
759            XTbl    -   array[0..M-1], values of X
760            YTbl    -   array[0..M-1,0..N-1], values of Y in X[i]
761            Rep     -   solver report:
762                        * Rep.TerminationType completetion code:
763                            * -2    X is not ordered  by  ascending/descending  or
764                                    there are non-distinct X[],  i.e.  X[i]=X[i+1]
765                            * -1    incorrect parameters were specified
766                            *  1    task has been solved
767                        * Rep.NFEV contains number of function calculations
768
769          -- ALGLIB --
770             Copyright 01.09.2009 by Bochkanov Sergey
771        *************************************************************************/
772        public static void odesolverresults(odesolverstate state,
773            ref int m,
774            ref double[] xtbl,
775            ref double[,] ytbl,
776            odesolverreport rep)
777        {
778            double v = 0;
779            int i = 0;
780            int i_ = 0;
781
782            m = 0;
783            xtbl = new double[0];
784            ytbl = new double[0,0];
785
786            rep.terminationtype = state.repterminationtype;
787            if( rep.terminationtype>0 )
788            {
789                m = state.m;
790                rep.nfev = state.repnfev;
791                xtbl = new double[state.m];
792                v = state.xscale;
793                for(i_=0; i_<=state.m-1;i_++)
794                {
795                    xtbl[i_] = v*state.xg[i_];
796                }
797                ytbl = new double[state.m, state.n];
798                for(i=0; i<=state.m-1; i++)
799                {
800                    for(i_=0; i_<=state.n-1;i_++)
801                    {
802                        ytbl[i,i_] = state.ytbl[i,i_];
803                    }
804                }
805            }
806            else
807            {
808                rep.nfev = 0;
809            }
810        }
811
812
813        /*************************************************************************
814        Internal initialization subroutine
815        *************************************************************************/
816        private static void odesolverinit(int solvertype,
817            double[] y,
818            int n,
819            double[] x,
820            int m,
821            double eps,
822            double h,
823            odesolverstate state)
824        {
825            int i = 0;
826            double v = 0;
827            int i_ = 0;
828
829           
830            //
831            // Prepare RComm
832            //
833            state.rstate.ia = new int[5+1];
834            state.rstate.ba = new bool[0+1];
835            state.rstate.ra = new double[5+1];
836            state.rstate.stage = -1;
837            state.needdy = false;
838           
839            //
840            // check parameters.
841            //
842            if( (n<=0 || m<1) || (double)(eps)==(double)(0) )
843            {
844                state.repterminationtype = -1;
845                return;
846            }
847            if( (double)(h)<(double)(0) )
848            {
849                h = -h;
850            }
851           
852            //
853            // quick exit if necessary.
854            // after this block we assume that M>1
855            //
856            if( m==1 )
857            {
858                state.repnfev = 0;
859                state.repterminationtype = 1;
860                state.ytbl = new double[1, n];
861                for(i_=0; i_<=n-1;i_++)
862                {
863                    state.ytbl[0,i_] = y[i_];
864                }
865                state.xg = new double[m];
866                for(i_=0; i_<=m-1;i_++)
867                {
868                    state.xg[i_] = x[i_];
869                }
870                return;
871            }
872           
873            //
874            // check again: correct order of X[]
875            //
876            if( (double)(x[1])==(double)(x[0]) )
877            {
878                state.repterminationtype = -2;
879                return;
880            }
881            for(i=1; i<=m-1; i++)
882            {
883                if( ((double)(x[1])>(double)(x[0]) && (double)(x[i])<=(double)(x[i-1])) || ((double)(x[1])<(double)(x[0]) && (double)(x[i])>=(double)(x[i-1])) )
884                {
885                    state.repterminationtype = -2;
886                    return;
887                }
888            }
889           
890            //
891            // auto-select H if necessary
892            //
893            if( (double)(h)==(double)(0) )
894            {
895                v = Math.Abs(x[1]-x[0]);
896                for(i=2; i<=m-1; i++)
897                {
898                    v = Math.Min(v, Math.Abs(x[i]-x[i-1]));
899                }
900                h = 0.001*v;
901            }
902           
903            //
904            // store parameters
905            //
906            state.n = n;
907            state.m = m;
908            state.h = h;
909            state.eps = Math.Abs(eps);
910            state.fraceps = (double)(eps)<(double)(0);
911            state.xg = new double[m];
912            for(i_=0; i_<=m-1;i_++)
913            {
914                state.xg[i_] = x[i_];
915            }
916            if( (double)(x[1])>(double)(x[0]) )
917            {
918                state.xscale = 1;
919            }
920            else
921            {
922                state.xscale = -1;
923                for(i_=0; i_<=m-1;i_++)
924                {
925                    state.xg[i_] = -1*state.xg[i_];
926                }
927            }
928            state.yc = new double[n];
929            for(i_=0; i_<=n-1;i_++)
930            {
931                state.yc[i_] = y[i_];
932            }
933            state.solvertype = solvertype;
934            state.repterminationtype = 0;
935           
936            //
937            // Allocate arrays
938            //
939            state.y = new double[n];
940            state.dy = new double[n];
941        }
942
943
944    }
945}
946
Note: See TracBrowser for help on using the repository browser.