[2806] | 1 | /*************************************************************************
|
---|
| 2 | Copyright 2009 by Sergey Bochkanov (ALGLIB project).
|
---|
| 3 |
|
---|
| 4 | >>> SOURCE LICENSE >>>
|
---|
| 5 | This program is free software; you can redistribute it and/or modify
|
---|
| 6 | it under the terms of the GNU General Public License as published by
|
---|
| 7 | the Free Software Foundation (www.fsf.org); either version 2 of the
|
---|
| 8 | License, or (at your option) any later version.
|
---|
| 9 |
|
---|
| 10 | This program is distributed in the hope that it will be useful,
|
---|
| 11 | but WITHOUT ANY WARRANTY; without even the implied warranty of
|
---|
| 12 | MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
---|
| 13 | GNU General Public License for more details.
|
---|
| 14 |
|
---|
| 15 | A copy of the GNU General Public License is available at
|
---|
| 16 | http://www.fsf.org/licensing/licenses
|
---|
| 17 |
|
---|
| 18 | >>> END OF LICENSE >>>
|
---|
| 19 | *************************************************************************/
|
---|
| 20 |
|
---|
| 21 | using System;
|
---|
| 22 |
|
---|
| 23 | namespace alglib
|
---|
| 24 | {
|
---|
| 25 | public class odesolver
|
---|
| 26 | {
|
---|
| 27 | public struct odesolverstate
|
---|
| 28 | {
|
---|
| 29 | public int n;
|
---|
| 30 | public int m;
|
---|
| 31 | public double xscale;
|
---|
| 32 | public double h;
|
---|
| 33 | public double eps;
|
---|
| 34 | public bool fraceps;
|
---|
| 35 | public double[] yc;
|
---|
| 36 | public double[] escale;
|
---|
| 37 | public double[] xg;
|
---|
| 38 | public int solvertype;
|
---|
| 39 | public double x;
|
---|
| 40 | public double[] y;
|
---|
| 41 | public double[] dy;
|
---|
| 42 | public double[,] ytbl;
|
---|
| 43 | public int repterminationtype;
|
---|
| 44 | public int repnfev;
|
---|
| 45 | public double[] yn;
|
---|
| 46 | public double[] yns;
|
---|
| 47 | public double[] rka;
|
---|
| 48 | public double[] rkc;
|
---|
| 49 | public double[] rkcs;
|
---|
| 50 | public double[,] rkb;
|
---|
| 51 | public double[,] rkk;
|
---|
| 52 | public AP.rcommstate rstate;
|
---|
| 53 | };
|
---|
| 54 |
|
---|
| 55 |
|
---|
| 56 | public struct odesolverreport
|
---|
| 57 | {
|
---|
| 58 | public int nfev;
|
---|
| 59 | public int terminationtype;
|
---|
| 60 | };
|
---|
| 61 |
|
---|
| 62 |
|
---|
| 63 |
|
---|
| 64 |
|
---|
| 65 | public const double odesolvermaxgrow = 3.0;
|
---|
| 66 | public const double odesolvermaxshrink = 10.0;
|
---|
| 67 |
|
---|
| 68 |
|
---|
| 69 | /*************************************************************************
|
---|
| 70 | Cash-Karp adaptive ODE solver.
|
---|
| 71 |
|
---|
| 72 | This subroutine solves ODE Y'=f(Y,x) with initial conditions Y(xs)=Ys
|
---|
| 73 | (here Y may be single variable or vector of N variables).
|
---|
| 74 |
|
---|
| 75 | INPUT PARAMETERS:
|
---|
| 76 | Y - initial conditions, array[0..N-1].
|
---|
| 77 | contains values of Y[] at X[0]
|
---|
| 78 | N - system size
|
---|
| 79 | X - points at which Y should be tabulated, array[0..M-1]
|
---|
| 80 | integrations starts at X[0], ends at X[M-1], intermediate
|
---|
| 81 | values at X[i] are returned too.
|
---|
| 82 | SHOULD BE ORDERED BY ASCENDING OR BY DESCENDING!!!!
|
---|
| 83 | M - number of intermediate points + first point + last point:
|
---|
| 84 | * M>2 means that you need both Y(X[M-1]) and M-2 values at
|
---|
| 85 | intermediate points
|
---|
| 86 | * M=2 means that you want just to integrate from X[0] to
|
---|
| 87 | X[1] and don't interested in intermediate values.
|
---|
| 88 | * M=1 means that you don't want to integrate :)
|
---|
| 89 | it is degenerate case, but it will be handled correctly.
|
---|
| 90 | * M<1 means error
|
---|
| 91 | Eps - tolerance (absolute/relative error on each step will be
|
---|
| 92 | less than Eps). When passing:
|
---|
| 93 | * Eps>0, it means desired ABSOLUTE error
|
---|
| 94 | * Eps<0, it means desired RELATIVE error. Relative errors
|
---|
| 95 | are calculated with respect to maximum values of Y seen
|
---|
| 96 | so far. Be careful to use this criterion when starting
|
---|
| 97 | from Y[] that are close to zero.
|
---|
| 98 | H - initial step lenth, it will be adjusted automatically
|
---|
| 99 | after the first step. If H=0, step will be selected
|
---|
| 100 | automatically (usualy it will be equal to 0.001 of
|
---|
| 101 | min(x[i]-x[j])).
|
---|
| 102 |
|
---|
| 103 | OUTPUT PARAMETERS
|
---|
| 104 | State - structure which stores algorithm state between subsequent
|
---|
| 105 | calls of OdeSolverIteration. Used for reverse communication.
|
---|
| 106 | This structure should be passed to the OdeSolverIteration
|
---|
| 107 | subroutine.
|
---|
| 108 |
|
---|
| 109 | SEE ALSO
|
---|
| 110 | AutoGKSmoothW, AutoGKSingular, AutoGKIteration, AutoGKResults.
|
---|
| 111 |
|
---|
| 112 |
|
---|
| 113 | -- ALGLIB --
|
---|
| 114 | Copyright 01.09.2009 by Bochkanov Sergey
|
---|
| 115 | *************************************************************************/
|
---|
| 116 | public static void odesolverrkck(ref double[] y,
|
---|
| 117 | int n,
|
---|
| 118 | ref double[] x,
|
---|
| 119 | int m,
|
---|
| 120 | double eps,
|
---|
| 121 | double h,
|
---|
| 122 | ref odesolverstate state)
|
---|
| 123 | {
|
---|
| 124 | odesolverinit(0, ref y, n, ref x, m, eps, h, ref state);
|
---|
| 125 | }
|
---|
| 126 |
|
---|
| 127 |
|
---|
| 128 | /*************************************************************************
|
---|
| 129 | One iteration of ODE solver.
|
---|
| 130 |
|
---|
| 131 | Called after inialization of State structure with OdeSolverXXX subroutine.
|
---|
| 132 | See HTML docs for examples.
|
---|
| 133 |
|
---|
| 134 | INPUT PARAMETERS:
|
---|
| 135 | State - structure which stores algorithm state between subsequent
|
---|
| 136 | calls and which is used for reverse communication. Must be
|
---|
| 137 | initialized with OdeSolverXXX() call first.
|
---|
| 138 |
|
---|
| 139 | If subroutine returned False, algorithm have finished its work.
|
---|
| 140 | If subroutine returned True, then user should:
|
---|
| 141 | * calculate F(State.X, State.Y)
|
---|
| 142 | * store it in State.DY
|
---|
| 143 | Here State.X is real, State.Y and State.DY are arrays[0..N-1] of reals.
|
---|
| 144 |
|
---|
| 145 | -- ALGLIB --
|
---|
| 146 | Copyright 01.09.2009 by Bochkanov Sergey
|
---|
| 147 | *************************************************************************/
|
---|
| 148 | public static bool odesolveriteration(ref odesolverstate state)
|
---|
| 149 | {
|
---|
| 150 | bool result = new bool();
|
---|
| 151 | int n = 0;
|
---|
| 152 | int m = 0;
|
---|
| 153 | int i = 0;
|
---|
| 154 | int j = 0;
|
---|
| 155 | int k = 0;
|
---|
| 156 | double xc = 0;
|
---|
| 157 | double v = 0;
|
---|
| 158 | double h = 0;
|
---|
| 159 | double h2 = 0;
|
---|
| 160 | bool gridpoint = new bool();
|
---|
| 161 | double err = 0;
|
---|
| 162 | double maxgrowpow = 0;
|
---|
| 163 | int klimit = 0;
|
---|
| 164 | int i_ = 0;
|
---|
| 165 |
|
---|
| 166 |
|
---|
| 167 | //
|
---|
| 168 | // Reverse communication preparations
|
---|
| 169 | // I know it looks ugly, but it works the same way
|
---|
| 170 | // anywhere from C++ to Python.
|
---|
| 171 | //
|
---|
| 172 | // This code initializes locals by:
|
---|
| 173 | // * random values determined during code
|
---|
| 174 | // generation - on first subroutine call
|
---|
| 175 | // * values from previous call - on subsequent calls
|
---|
| 176 | //
|
---|
| 177 | if( state.rstate.stage>=0 )
|
---|
| 178 | {
|
---|
| 179 | n = state.rstate.ia[0];
|
---|
| 180 | m = state.rstate.ia[1];
|
---|
| 181 | i = state.rstate.ia[2];
|
---|
| 182 | j = state.rstate.ia[3];
|
---|
| 183 | k = state.rstate.ia[4];
|
---|
| 184 | klimit = state.rstate.ia[5];
|
---|
| 185 | gridpoint = state.rstate.ba[0];
|
---|
| 186 | xc = state.rstate.ra[0];
|
---|
| 187 | v = state.rstate.ra[1];
|
---|
| 188 | h = state.rstate.ra[2];
|
---|
| 189 | h2 = state.rstate.ra[3];
|
---|
| 190 | err = state.rstate.ra[4];
|
---|
| 191 | maxgrowpow = state.rstate.ra[5];
|
---|
| 192 | }
|
---|
| 193 | else
|
---|
| 194 | {
|
---|
| 195 | n = -983;
|
---|
| 196 | m = -989;
|
---|
| 197 | i = -834;
|
---|
| 198 | j = 900;
|
---|
| 199 | k = -287;
|
---|
| 200 | klimit = 364;
|
---|
| 201 | gridpoint = false;
|
---|
| 202 | xc = -338;
|
---|
| 203 | v = -686;
|
---|
| 204 | h = 912;
|
---|
| 205 | h2 = 585;
|
---|
| 206 | err = 497;
|
---|
| 207 | maxgrowpow = -271;
|
---|
| 208 | }
|
---|
| 209 | if( state.rstate.stage==0 )
|
---|
| 210 | {
|
---|
| 211 | goto lbl_0;
|
---|
| 212 | }
|
---|
| 213 |
|
---|
| 214 | //
|
---|
| 215 | // Routine body
|
---|
| 216 | //
|
---|
| 217 |
|
---|
| 218 | //
|
---|
| 219 | // prepare
|
---|
| 220 | //
|
---|
| 221 | if( state.repterminationtype!=0 )
|
---|
| 222 | {
|
---|
| 223 | result = false;
|
---|
| 224 | return result;
|
---|
| 225 | }
|
---|
| 226 | n = state.n;
|
---|
| 227 | m = state.m;
|
---|
| 228 | h = state.h;
|
---|
| 229 | state.y = new double[n];
|
---|
| 230 | state.dy = new double[n];
|
---|
| 231 | maxgrowpow = Math.Pow(odesolvermaxgrow, 5);
|
---|
| 232 | state.repnfev = 0;
|
---|
| 233 |
|
---|
| 234 | //
|
---|
| 235 | // some preliminary checks for internal errors
|
---|
| 236 | // after this we assume that H>0 and M>1
|
---|
| 237 | //
|
---|
| 238 | System.Diagnostics.Debug.Assert((double)(state.h)>(double)(0), "ODESolver: internal error");
|
---|
| 239 | System.Diagnostics.Debug.Assert(m>1, "ODESolverIteration: internal error");
|
---|
| 240 |
|
---|
| 241 | //
|
---|
| 242 | // choose solver
|
---|
| 243 | //
|
---|
| 244 | if( state.solvertype!=0 )
|
---|
| 245 | {
|
---|
| 246 | goto lbl_1;
|
---|
| 247 | }
|
---|
| 248 |
|
---|
| 249 | //
|
---|
| 250 | // Cask-Karp solver
|
---|
| 251 | // Prepare coefficients table.
|
---|
| 252 | // Check it for errors
|
---|
| 253 | //
|
---|
| 254 | state.rka = new double[6];
|
---|
| 255 | state.rka[0] = 0;
|
---|
| 256 | state.rka[1] = (double)(1)/(double)(5);
|
---|
| 257 | state.rka[2] = (double)(3)/(double)(10);
|
---|
| 258 | state.rka[3] = (double)(3)/(double)(5);
|
---|
| 259 | state.rka[4] = 1;
|
---|
| 260 | state.rka[5] = (double)(7)/(double)(8);
|
---|
| 261 | state.rkb = new double[6, 5];
|
---|
| 262 | state.rkb[1,0] = (double)(1)/(double)(5);
|
---|
| 263 | state.rkb[2,0] = (double)(3)/(double)(40);
|
---|
| 264 | state.rkb[2,1] = (double)(9)/(double)(40);
|
---|
| 265 | state.rkb[3,0] = (double)(3)/(double)(10);
|
---|
| 266 | state.rkb[3,1] = -((double)(9)/(double)(10));
|
---|
| 267 | state.rkb[3,2] = (double)(6)/(double)(5);
|
---|
| 268 | state.rkb[4,0] = -((double)(11)/(double)(54));
|
---|
| 269 | state.rkb[4,1] = (double)(5)/(double)(2);
|
---|
| 270 | state.rkb[4,2] = -((double)(70)/(double)(27));
|
---|
| 271 | state.rkb[4,3] = (double)(35)/(double)(27);
|
---|
| 272 | state.rkb[5,0] = (double)(1631)/(double)(55296);
|
---|
| 273 | state.rkb[5,1] = (double)(175)/(double)(512);
|
---|
| 274 | state.rkb[5,2] = (double)(575)/(double)(13824);
|
---|
| 275 | state.rkb[5,3] = (double)(44275)/(double)(110592);
|
---|
| 276 | state.rkb[5,4] = (double)(253)/(double)(4096);
|
---|
| 277 | state.rkc = new double[6];
|
---|
| 278 | state.rkc[0] = (double)(37)/(double)(378);
|
---|
| 279 | state.rkc[1] = 0;
|
---|
| 280 | state.rkc[2] = (double)(250)/(double)(621);
|
---|
| 281 | state.rkc[3] = (double)(125)/(double)(594);
|
---|
| 282 | state.rkc[4] = 0;
|
---|
| 283 | state.rkc[5] = (double)(512)/(double)(1771);
|
---|
| 284 | state.rkcs = new double[6];
|
---|
| 285 | state.rkcs[0] = (double)(2825)/(double)(27648);
|
---|
| 286 | state.rkcs[1] = 0;
|
---|
| 287 | state.rkcs[2] = (double)(18575)/(double)(48384);
|
---|
| 288 | state.rkcs[3] = (double)(13525)/(double)(55296);
|
---|
| 289 | state.rkcs[4] = (double)(277)/(double)(14336);
|
---|
| 290 | state.rkcs[5] = (double)(1)/(double)(4);
|
---|
| 291 | state.rkk = new double[6, n];
|
---|
| 292 |
|
---|
| 293 | //
|
---|
| 294 | // Main cycle consists of two iterations:
|
---|
| 295 | // * outer where we travel from X[i-1] to X[i]
|
---|
| 296 | // * inner where we travel inside [X[i-1],X[i]]
|
---|
| 297 | //
|
---|
| 298 | state.ytbl = new double[m, n];
|
---|
| 299 | state.escale = new double[n];
|
---|
| 300 | state.yn = new double[n];
|
---|
| 301 | state.yns = new double[n];
|
---|
| 302 | xc = state.xg[0];
|
---|
| 303 | for(i_=0; i_<=n-1;i_++)
|
---|
| 304 | {
|
---|
| 305 | state.ytbl[0,i_] = state.yc[i_];
|
---|
| 306 | }
|
---|
| 307 | for(j=0; j<=n-1; j++)
|
---|
| 308 | {
|
---|
| 309 | state.escale[j] = 0;
|
---|
| 310 | }
|
---|
| 311 | i = 1;
|
---|
| 312 | lbl_3:
|
---|
| 313 | if( i>m-1 )
|
---|
| 314 | {
|
---|
| 315 | goto lbl_5;
|
---|
| 316 | }
|
---|
| 317 |
|
---|
| 318 | //
|
---|
| 319 | // begin inner iteration
|
---|
| 320 | //
|
---|
| 321 | lbl_6:
|
---|
| 322 | if( false )
|
---|
| 323 | {
|
---|
| 324 | goto lbl_7;
|
---|
| 325 | }
|
---|
| 326 |
|
---|
| 327 | //
|
---|
| 328 | // truncate step if needed (beyond right boundary).
|
---|
| 329 | // determine should we store X or not
|
---|
| 330 | //
|
---|
| 331 | if( (double)(xc+h)>=(double)(state.xg[i]) )
|
---|
| 332 | {
|
---|
| 333 | h = state.xg[i]-xc;
|
---|
| 334 | gridpoint = true;
|
---|
| 335 | }
|
---|
| 336 | else
|
---|
| 337 | {
|
---|
| 338 | gridpoint = false;
|
---|
| 339 | }
|
---|
| 340 |
|
---|
| 341 | //
|
---|
| 342 | // Update error scale maximums
|
---|
| 343 | //
|
---|
| 344 | // These maximums are initialized by zeros,
|
---|
| 345 | // then updated every iterations.
|
---|
| 346 | //
|
---|
| 347 | for(j=0; j<=n-1; j++)
|
---|
| 348 | {
|
---|
| 349 | state.escale[j] = Math.Max(state.escale[j], Math.Abs(state.yc[j]));
|
---|
| 350 | }
|
---|
| 351 |
|
---|
| 352 | //
|
---|
| 353 | // make one step:
|
---|
| 354 | // 1. calculate all info needed to do step
|
---|
| 355 | // 2. update errors scale maximums using values/derivatives
|
---|
| 356 | // obtained during (1)
|
---|
| 357 | //
|
---|
| 358 | // Take into account that we use scaling of X to reduce task
|
---|
| 359 | // to the form where x[0] < x[1] < ... < x[n-1]. So X is
|
---|
| 360 | // replaced by x=xscale*t, and dy/dx=f(y,x) is replaced
|
---|
| 361 | // by dy/dt=xscale*f(y,xscale*t).
|
---|
| 362 | //
|
---|
| 363 | for(i_=0; i_<=n-1;i_++)
|
---|
| 364 | {
|
---|
| 365 | state.yn[i_] = state.yc[i_];
|
---|
| 366 | }
|
---|
| 367 | for(i_=0; i_<=n-1;i_++)
|
---|
| 368 | {
|
---|
| 369 | state.yns[i_] = state.yc[i_];
|
---|
| 370 | }
|
---|
| 371 | k = 0;
|
---|
| 372 | lbl_8:
|
---|
| 373 | if( k>5 )
|
---|
| 374 | {
|
---|
| 375 | goto lbl_10;
|
---|
| 376 | }
|
---|
| 377 |
|
---|
| 378 | //
|
---|
| 379 | // prepare data for the next update of YN/YNS
|
---|
| 380 | //
|
---|
| 381 | state.x = state.xscale*(xc+state.rka[k]*h);
|
---|
| 382 | for(i_=0; i_<=n-1;i_++)
|
---|
| 383 | {
|
---|
| 384 | state.y[i_] = state.yc[i_];
|
---|
| 385 | }
|
---|
| 386 | for(j=0; j<=k-1; j++)
|
---|
| 387 | {
|
---|
| 388 | v = state.rkb[k,j];
|
---|
| 389 | for(i_=0; i_<=n-1;i_++)
|
---|
| 390 | {
|
---|
| 391 | state.y[i_] = state.y[i_] + v*state.rkk[j,i_];
|
---|
| 392 | }
|
---|
| 393 | }
|
---|
| 394 | state.rstate.stage = 0;
|
---|
| 395 | goto lbl_rcomm;
|
---|
| 396 | lbl_0:
|
---|
| 397 | state.repnfev = state.repnfev+1;
|
---|
| 398 | v = h*state.xscale;
|
---|
| 399 | for(i_=0; i_<=n-1;i_++)
|
---|
| 400 | {
|
---|
| 401 | state.rkk[k,i_] = v*state.dy[i_];
|
---|
| 402 | }
|
---|
| 403 |
|
---|
| 404 | //
|
---|
| 405 | // update YN/YNS
|
---|
| 406 | //
|
---|
| 407 | v = state.rkc[k];
|
---|
| 408 | for(i_=0; i_<=n-1;i_++)
|
---|
| 409 | {
|
---|
| 410 | state.yn[i_] = state.yn[i_] + v*state.rkk[k,i_];
|
---|
| 411 | }
|
---|
| 412 | v = state.rkcs[k];
|
---|
| 413 | for(i_=0; i_<=n-1;i_++)
|
---|
| 414 | {
|
---|
| 415 | state.yns[i_] = state.yns[i_] + v*state.rkk[k,i_];
|
---|
| 416 | }
|
---|
| 417 | k = k+1;
|
---|
| 418 | goto lbl_8;
|
---|
| 419 | lbl_10:
|
---|
| 420 |
|
---|
| 421 | //
|
---|
| 422 | // estimate error
|
---|
| 423 | //
|
---|
| 424 | err = 0;
|
---|
| 425 | for(j=0; j<=n-1; j++)
|
---|
| 426 | {
|
---|
| 427 | if( !state.fraceps )
|
---|
| 428 | {
|
---|
| 429 |
|
---|
| 430 | //
|
---|
| 431 | // absolute error is estimated
|
---|
| 432 | //
|
---|
| 433 | err = Math.Max(err, Math.Abs(state.yn[j]-state.yns[j]));
|
---|
| 434 | }
|
---|
| 435 | else
|
---|
| 436 | {
|
---|
| 437 |
|
---|
| 438 | //
|
---|
| 439 | // Relative error is estimated
|
---|
| 440 | //
|
---|
| 441 | v = state.escale[j];
|
---|
| 442 | if( (double)(v)==(double)(0) )
|
---|
| 443 | {
|
---|
| 444 | v = 1;
|
---|
| 445 | }
|
---|
| 446 | err = Math.Max(err, Math.Abs(state.yn[j]-state.yns[j])/v);
|
---|
| 447 | }
|
---|
| 448 | }
|
---|
| 449 |
|
---|
| 450 | //
|
---|
| 451 | // calculate new step, restart if necessary
|
---|
| 452 | //
|
---|
| 453 | if( (double)(maxgrowpow*err)<=(double)(state.eps) )
|
---|
| 454 | {
|
---|
| 455 | h2 = odesolvermaxgrow*h;
|
---|
| 456 | }
|
---|
| 457 | else
|
---|
| 458 | {
|
---|
| 459 | h2 = h*Math.Pow(state.eps/err, 0.2);
|
---|
| 460 | }
|
---|
| 461 | if( (double)(h2)<(double)(h/odesolvermaxshrink) )
|
---|
| 462 | {
|
---|
| 463 | h2 = h/odesolvermaxshrink;
|
---|
| 464 | }
|
---|
| 465 | if( (double)(err)>(double)(state.eps) )
|
---|
| 466 | {
|
---|
| 467 | h = h2;
|
---|
| 468 | goto lbl_6;
|
---|
| 469 | }
|
---|
| 470 |
|
---|
| 471 | //
|
---|
| 472 | // advance position
|
---|
| 473 | //
|
---|
| 474 | xc = xc+h;
|
---|
| 475 | for(i_=0; i_<=n-1;i_++)
|
---|
| 476 | {
|
---|
| 477 | state.yc[i_] = state.yn[i_];
|
---|
| 478 | }
|
---|
| 479 |
|
---|
| 480 | //
|
---|
| 481 | // update H
|
---|
| 482 | //
|
---|
| 483 | h = h2;
|
---|
| 484 |
|
---|
| 485 | //
|
---|
| 486 | // break on grid point
|
---|
| 487 | //
|
---|
| 488 | if( gridpoint )
|
---|
| 489 | {
|
---|
| 490 | goto lbl_7;
|
---|
| 491 | }
|
---|
| 492 | goto lbl_6;
|
---|
| 493 | lbl_7:
|
---|
| 494 |
|
---|
| 495 | //
|
---|
| 496 | // save result
|
---|
| 497 | //
|
---|
| 498 | for(i_=0; i_<=n-1;i_++)
|
---|
| 499 | {
|
---|
| 500 | state.ytbl[i,i_] = state.yc[i_];
|
---|
| 501 | }
|
---|
| 502 | i = i+1;
|
---|
| 503 | goto lbl_3;
|
---|
| 504 | lbl_5:
|
---|
| 505 | state.repterminationtype = 1;
|
---|
| 506 | result = false;
|
---|
| 507 | return result;
|
---|
| 508 | lbl_1:
|
---|
| 509 | result = false;
|
---|
| 510 | return result;
|
---|
| 511 |
|
---|
| 512 | //
|
---|
| 513 | // Saving state
|
---|
| 514 | //
|
---|
| 515 | lbl_rcomm:
|
---|
| 516 | result = true;
|
---|
| 517 | state.rstate.ia[0] = n;
|
---|
| 518 | state.rstate.ia[1] = m;
|
---|
| 519 | state.rstate.ia[2] = i;
|
---|
| 520 | state.rstate.ia[3] = j;
|
---|
| 521 | state.rstate.ia[4] = k;
|
---|
| 522 | state.rstate.ia[5] = klimit;
|
---|
| 523 | state.rstate.ba[0] = gridpoint;
|
---|
| 524 | state.rstate.ra[0] = xc;
|
---|
| 525 | state.rstate.ra[1] = v;
|
---|
| 526 | state.rstate.ra[2] = h;
|
---|
| 527 | state.rstate.ra[3] = h2;
|
---|
| 528 | state.rstate.ra[4] = err;
|
---|
| 529 | state.rstate.ra[5] = maxgrowpow;
|
---|
| 530 | return result;
|
---|
| 531 | }
|
---|
| 532 |
|
---|
| 533 |
|
---|
| 534 | /*************************************************************************
|
---|
| 535 | ODE solver results
|
---|
| 536 |
|
---|
| 537 | Called after OdeSolverIteration returned False.
|
---|
| 538 |
|
---|
| 539 | INPUT PARAMETERS:
|
---|
| 540 | State - algorithm state (used by OdeSolverIteration).
|
---|
| 541 |
|
---|
| 542 | OUTPUT PARAMETERS:
|
---|
| 543 | M - number of tabulated values, M>=1
|
---|
| 544 | XTbl - array[0..M-1], values of X
|
---|
| 545 | YTbl - array[0..M-1,0..N-1], values of Y in X[i]
|
---|
| 546 | Rep - solver report:
|
---|
| 547 | * Rep.TerminationType completetion code:
|
---|
| 548 | * -2 X is not ordered by ascending/descending or
|
---|
| 549 | there are non-distinct X[], i.e. X[i]=X[i+1]
|
---|
| 550 | * -1 incorrect parameters were specified
|
---|
| 551 | * 1 task has been solved
|
---|
| 552 | * Rep.NFEV contains number of function calculations
|
---|
| 553 |
|
---|
| 554 | -- ALGLIB --
|
---|
| 555 | Copyright 01.09.2009 by Bochkanov Sergey
|
---|
| 556 | *************************************************************************/
|
---|
| 557 | public static void odesolverresults(ref odesolverstate state,
|
---|
| 558 | ref int m,
|
---|
| 559 | ref double[] xtbl,
|
---|
| 560 | ref double[,] ytbl,
|
---|
| 561 | ref odesolverreport rep)
|
---|
| 562 | {
|
---|
| 563 | double v = 0;
|
---|
| 564 | int i = 0;
|
---|
| 565 | int i_ = 0;
|
---|
| 566 |
|
---|
| 567 | rep.terminationtype = state.repterminationtype;
|
---|
| 568 | if( rep.terminationtype>0 )
|
---|
| 569 | {
|
---|
| 570 | m = state.m;
|
---|
| 571 | rep.nfev = state.repnfev;
|
---|
| 572 | xtbl = new double[state.m];
|
---|
| 573 | v = state.xscale;
|
---|
| 574 | for(i_=0; i_<=state.m-1;i_++)
|
---|
| 575 | {
|
---|
| 576 | xtbl[i_] = v*state.xg[i_];
|
---|
| 577 | }
|
---|
| 578 | ytbl = new double[state.m, state.n];
|
---|
| 579 | for(i=0; i<=state.m-1; i++)
|
---|
| 580 | {
|
---|
| 581 | for(i_=0; i_<=state.n-1;i_++)
|
---|
| 582 | {
|
---|
| 583 | ytbl[i,i_] = state.ytbl[i,i_];
|
---|
| 584 | }
|
---|
| 585 | }
|
---|
| 586 | }
|
---|
| 587 | else
|
---|
| 588 | {
|
---|
| 589 | rep.nfev = 0;
|
---|
| 590 | }
|
---|
| 591 | }
|
---|
| 592 |
|
---|
| 593 |
|
---|
| 594 | /*************************************************************************
|
---|
| 595 | Internal initialization subroutine
|
---|
| 596 | *************************************************************************/
|
---|
| 597 | private static void odesolverinit(int solvertype,
|
---|
| 598 | ref double[] y,
|
---|
| 599 | int n,
|
---|
| 600 | ref double[] x,
|
---|
| 601 | int m,
|
---|
| 602 | double eps,
|
---|
| 603 | double h,
|
---|
| 604 | ref odesolverstate state)
|
---|
| 605 | {
|
---|
| 606 | int i = 0;
|
---|
| 607 | double v = 0;
|
---|
| 608 | int i_ = 0;
|
---|
| 609 |
|
---|
| 610 |
|
---|
| 611 | //
|
---|
| 612 | // Prepare RComm
|
---|
| 613 | //
|
---|
| 614 | state.rstate.ia = new int[5+1];
|
---|
| 615 | state.rstate.ba = new bool[0+1];
|
---|
| 616 | state.rstate.ra = new double[5+1];
|
---|
| 617 | state.rstate.stage = -1;
|
---|
| 618 |
|
---|
| 619 | //
|
---|
| 620 | // check parameters.
|
---|
| 621 | //
|
---|
| 622 | if( n<=0 | m<1 | (double)(eps)==(double)(0) )
|
---|
| 623 | {
|
---|
| 624 | state.repterminationtype = -1;
|
---|
| 625 | return;
|
---|
| 626 | }
|
---|
| 627 | if( (double)(h)<(double)(0) )
|
---|
| 628 | {
|
---|
| 629 | h = -h;
|
---|
| 630 | }
|
---|
| 631 |
|
---|
| 632 | //
|
---|
| 633 | // quick exit if necessary.
|
---|
| 634 | // after this block we assume that M>1
|
---|
| 635 | //
|
---|
| 636 | if( m==1 )
|
---|
| 637 | {
|
---|
| 638 | state.repnfev = 0;
|
---|
| 639 | state.repterminationtype = 1;
|
---|
| 640 | state.ytbl = new double[1, n];
|
---|
| 641 | for(i_=0; i_<=n-1;i_++)
|
---|
| 642 | {
|
---|
| 643 | state.ytbl[0,i_] = y[i_];
|
---|
| 644 | }
|
---|
| 645 | state.xg = new double[m];
|
---|
| 646 | for(i_=0; i_<=m-1;i_++)
|
---|
| 647 | {
|
---|
| 648 | state.xg[i_] = x[i_];
|
---|
| 649 | }
|
---|
| 650 | return;
|
---|
| 651 | }
|
---|
| 652 |
|
---|
| 653 | //
|
---|
| 654 | // check again: correct order of X[]
|
---|
| 655 | //
|
---|
| 656 | if( (double)(x[1])==(double)(x[0]) )
|
---|
| 657 | {
|
---|
| 658 | state.repterminationtype = -2;
|
---|
| 659 | return;
|
---|
| 660 | }
|
---|
| 661 | for(i=1; i<=m-1; i++)
|
---|
| 662 | {
|
---|
| 663 | 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]) )
|
---|
| 664 | {
|
---|
| 665 | state.repterminationtype = -2;
|
---|
| 666 | return;
|
---|
| 667 | }
|
---|
| 668 | }
|
---|
| 669 |
|
---|
| 670 | //
|
---|
| 671 | // auto-select H if necessary
|
---|
| 672 | //
|
---|
| 673 | if( (double)(h)==(double)(0) )
|
---|
| 674 | {
|
---|
| 675 | v = Math.Abs(x[1]-x[0]);
|
---|
| 676 | for(i=2; i<=m-1; i++)
|
---|
| 677 | {
|
---|
| 678 | v = Math.Min(v, Math.Abs(x[i]-x[i-1]));
|
---|
| 679 | }
|
---|
| 680 | h = 0.001*v;
|
---|
| 681 | }
|
---|
| 682 |
|
---|
| 683 | //
|
---|
| 684 | // store parameters
|
---|
| 685 | //
|
---|
| 686 | state.n = n;
|
---|
| 687 | state.m = m;
|
---|
| 688 | state.h = h;
|
---|
| 689 | state.eps = Math.Abs(eps);
|
---|
| 690 | state.fraceps = (double)(eps)<(double)(0);
|
---|
| 691 | state.xg = new double[m];
|
---|
| 692 | for(i_=0; i_<=m-1;i_++)
|
---|
| 693 | {
|
---|
| 694 | state.xg[i_] = x[i_];
|
---|
| 695 | }
|
---|
| 696 | if( (double)(x[1])>(double)(x[0]) )
|
---|
| 697 | {
|
---|
| 698 | state.xscale = 1;
|
---|
| 699 | }
|
---|
| 700 | else
|
---|
| 701 | {
|
---|
| 702 | state.xscale = -1;
|
---|
| 703 | for(i_=0; i_<=m-1;i_++)
|
---|
| 704 | {
|
---|
| 705 | state.xg[i_] = -1*state.xg[i_];
|
---|
| 706 | }
|
---|
| 707 | }
|
---|
| 708 | state.yc = new double[n];
|
---|
| 709 | for(i_=0; i_<=n-1;i_++)
|
---|
| 710 | {
|
---|
| 711 | state.yc[i_] = y[i_];
|
---|
| 712 | }
|
---|
| 713 | state.solvertype = solvertype;
|
---|
| 714 | state.repterminationtype = 0;
|
---|
| 715 | }
|
---|
| 716 | }
|
---|
| 717 | }
|
---|