/[ascend]/trunk/models/ivpNondimensional/ivpStepN.a4c
ViewVC logotype

Contents of /trunk/models/ivpNondimensional/ivpStepN.a4c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 5 - (show annotations) (download) (as text)
Sun Nov 7 19:45:14 2004 UTC (19 years, 7 months ago) by aw0a
File MIME type: text/x-ascend
File size: 45099 byte(s)
adding ivp models to model library
1 REQUIRE "atoms.a4l";
2
3 (*
4 * ivpStep.a4c
5 * by Arthur W. Westerberg
6 * Part of the ASCEND Library
7 * $Date: 02/07/28 $
8 * $Revision: 1.7 $
9 * $Author: mthomas $
10 * $Source: /afs/cs.cmu.edu/project/ascend/Repository/models/ivpsystem.a4l,v $
11 *
12 * This file is part of the ASCEND Modeling Library.
13 *
14 * Copyright (C) 1994 - 1998 Carnegie Mellon University
15 *
16 * The ASCEND Modeling Library is free software; you can redistribute
17 * it and/or modify it under the terms of the GNU General Public
18 * License as published by the Free Software Foundation; either
19 * version 2 of the License, or (at your option) any later version.
20 *
21 * The ASCEND Modeling Library is distributed in hope that it
22 * will be useful, but WITHOUT ANY WARRANTY; without even the implied
23 * warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
24 * See the GNU General Public License for more details.
25 *
26 * You should have received a copy of the GNU General Public License
27 * along with the program; if not, write to the Free Software
28 * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139 USA. Check
29 * the file named COPYING.
30 *)
31
32 (*============================================================================
33
34 I V P S T E P . A 4 C
35 -----------------------------
36
37 AUTHOR: Arthur W. Westerberg
38
39 DATES: 07/2002 - Original code (AWW)
40 07/2003 - Tested code. Works for simple test problem.
41 (AWW)
42 12/2003 - Added Tcl script that runs complete simulation
43 (AWW)
44
45 CONTENTS: Models for the numerical integration equations for
46 the two multistep methods: a bdf for stiff problems
47 and the Adams Moulton for non-stiff problems. There
48 is also a framework within which one creates the
49 model for the physical system. These models are for
50 taking one step of the independent variable.
51
52 ============================================================================*)
53
54
55 (* ---------------------------------------------------------- *)
56
57 MODEL ivpBase;
58
59 (* ----- ivpBase -----
60 The base model for solving DAE initial value problems using the
61 multistep methods: (a) a bdf method for stiff problems and (b)
62 the Adams Moulton method for non-stiff problems. A DAE ivp
63 model has two types of sets of equations:
64
65 (1) the model of the physical system at the k-th step
66
67 f(dydx(k), y(k), z(k), u(k), x(k)) = 0
68 h(y(k), z(k), u(k), x(k)) = 0
69 u(k) = u(x(k))
70
71 where y are the state, z the algebraic and u the forcing
72 variable vectors and x the scalar independent variable, and
73 where f is a vector of functions equal to the number of
74 states, h is a vector of algebraic equations equal in number
75 to the algebraic variables, and where one typically
76 specifies the forcing variables as functions only of the
77 independent variable.
78
79 (2) the numerical multistep integration equations
80
81 y(k) = y(k-1) + s(x[k]-x[k-1];
82 y[k-1],y[k-2],...;
83 dydx[k],dydx[k-1],...)
84
85 These equations relate the states at the current point to
86 the time derivatives at current and past points and states
87 at past points. They are not related to the model of the
88 physical system but only to the numerical method being used.
89
90 A common solving approach would be to solve the integration
91 equations in an outside loop, and, for each iteration of that
92 outside loop, use a Newton-based method to solve the physical
93 system model for the state derivatives and the algebraic
94 variables. Often this approach solves the outer loop using a
95 simple substitution method - which can be very slow to
96 converge. The following figure illustrates.
97
98
99 ------------------
100 u(x) ->---| (inner loop) |---- z -->
101 | solve model eqns |
102 ------>---| |-- dydx-->-
103 | ------------------ |
104 | |
105 | ------------------ |
106 | | (outer loop) | |
107 --- y,x =-| iterate inte- |---<-------
108 | gration eqns |
109 | |
110 ------------------
111
112
113 By creating both the model and integration equations through
114 ASCEND models, ASCEND can simultaneously converge the model
115 and integration equations using its sparse Newton-based solver
116 - thus making the solving much more efficient. We wrote this
117 model precisely to gain this efficiency.
118
119
120 The integration module - degrees of freedom considerations
121
122 For the q-th order bdf and Adams-Moulton (am) methods used
123 here, we write a q-th order polynomial in the independent
124 variable x for each state variable y[k][i]. This polynomial is
125 written in the form of a Taylor Series expanded around the
126 current point x[k].
127
128 A q-th order polynomial has the q+1 polynomial coefficients in
129 it a[k][i]/i!. The a[k][i] are estimate the q-th derivative of
130 the function at the current point, with the q-th term being an
131 estimate of the error in the fit. We determine these
132 coefficients by fitting this polynomial and/or the first
133 derivative of this polynomial to past y[j][i] and dy/dt[j][i],
134 j<k, data, respectively, and to the computed value for
135 dy/dt[k][i] at the current point. We would need to write
136 precisely q+1 such equations for this fitting. The bdf and am
137 methods differ in which equations we write to do this fitting.
138 See comments in later in the code for the actual polynomials we
139 write.
140
141 We then use this polynomial evaluated at the current point to
142 estimate the state y[k][i] at the current point x[k].
143
144 Looking back at the two boxes in the previous figure, the lower
145 integration box will have these q+1 equations in it, tabulated
146 values for previous points, and the means to use past data and
147 the current estimate for dy/dt[k][i] (the input we show above
148 for this box) to fit the coefficients for each state variable.
149 It will also write the polynomial one additional time to
150 provide the estimate for y[k][i] at the current point, which is
151 the output from this box. It thus contains q+2 polynomial
152 equations in it, in terms of q+1 unknown coefficients for each
153 state. The model equations introduce the states y[k][i] and their
154 derivatives dy/dx[k][i] for the current point and the equation to
155 compute the latter in terms of the former, often in the form
156
157 dy/dt[k][i] = f(all y[k],t).
158
159 So it introduces, in principle two variables and one equation
160 for each state.
161
162 To summarize, the integration box introduces q+1 equations and
163 q+2 variables per state while the model box introduces 1
164 equation and 2 new variables. The model is in square and
165 solvable.
166
167
168 Predictor behavior
169
170 The polynomial we fit for the current point can easily be used
171 to predict y[k+1][i], the values for the states i at the next
172 time step. We can also create polynomials to predict any of
173 the non state variables in the problem. We can include these
174 non state variables in the integration model and use the
175 current and q past values to determine the coefficients for a
176 Taylor Series polynomial that passes through these points.
177 This polynomial would only be useful for predicting as it play
178 no role in solving the problem. But often predicting one of
179 the algebraic variables might improve the chance to converge
180 the equations.
181
182
183 Error control
184
185 The package checks periodically to see discover the size step
186 it should take to keep the error at that the user has
187 prescribed. It can also change the order of the polynomial up
188 or down by one if such a move will allow for a larger step.
189
190 We write the polynomial in the form of a Taylor Series around
191 the current point. So the error in estimating the current
192 point, as a function of the step size we take, is the highest
193 order term in that series. We back compute the step size
194 needed to obtain the error we allow and use it as the step size
195 for the next integration step. We keep this step size for the
196 next q steps before checking again. Changing the order of the
197 polynomial can also affect the error. This package will change
198 the order up or down by one (i.e., to q-1 or q+1 (keeping in
199 the limits of 1 and 6) if that change allows it to take a
200 larger step. The next to last term in the Taylor Series
201 estimates the error when using a polynomial of order q-1. We
202 estimate the derivative for order q+1 as the difference in the
203 q-th derivative between the current point and the previous
204 point divided by the change in x between these two points -
205 i.e., a divided difference which is essentially a free
206 calculation.
207
208
209 Stopping conditions
210
211 We have created this package so one can stop whenever any of
212 the variables we are predicting (all the states and whichever
213 of the non states we select) crosses through zero. In taking a
214 normal step in this package, we use the error control criterion
215 to fix the step size in the independent variable, x, prior to
216 taking the next step. If after we solve the model for this
217 step, any one of the variables we have selected to use as
218 stopping conditions changes sign, we fix its value to zero and
219 resolve the entire model to compute the step size. As more
220 than one such variable may change sign in this final interval,
221 we test after each solution to see if another variable changed
222 sign in the reduced interval we just computed. If so, we
223 repeat by fixing this new variable to zero and resolving. In
224 this manner we will find the variable that changed sign at the
225 earliest point in the final interval.
226
227 *)
228
229 (* ----- ivpBase -----
230 TO USE THIS PACKAGE, the user must provide three models that
231 are refinements of models we provide here. The first two
232 will be refinements of integrationPoint, while the third must be
233 a refinement of ivpStep. See the comments for those models.
234 Also see the model ivpStep.testModel.a4c, which is a prototype
235 for what the user must write.
236
237 NOTE: THE MODELS IN THIS LIBRARY ARE WRITTEN USING DIMENSIONLESS
238 QUANTITIES. THE USER MUST CONVERT ANY DIMENSIONAL STATES
239 AND PREDICTED ALGEBRAIC VARIABLES, AND ALL STOPPING CONDITIONS,
240 ERRORS, ETC., TO DIMENSIONLESS VARIABLES FOR TYING INTO THESE
241 MODELS.
242
243 *)
244
245
246 END ivpBase;
247
248 (* ---------------------------------------------------------- *)
249
250 MODEL factorial(
251 N IS_A integer_constant;
252 )
253 WHERE(
254 N>=0;
255 N<=12; (* Note that 12! = 0.48E9, 2^31=2.1E9 *)
256 ) REFINES ivpBase;
257
258 (* create a vector of factorials, where fac[i] contains
259 i!, where i goes from 0 to N *)
260
261 fac[0..N] IS_A factor;
262
263 METHODS
264
265 (* ----- factorial ----- *)
266
267 METHOD default_self;
268 fac[0] := 1;
269
270 (* if N=0, the following loop produces no code *)
271 FOR i IN [1..N] DO
272 fac[i] := i*fac[i-1];
273 END FOR;
274 END default_self;
275
276 (* ----- factorial ----- *)
277
278 METHOD specify;
279 fac[0..N].fixed := TRUE;
280 END specify;
281
282 END factorial;
283
284 (* ---------------------------------------------------------- *)
285
286 MODEL pTrow(
287 set IS_A set OF integer_constant;
288 ) REFINES ivpBase;
289
290 (* This model is needed and instanced in
291 "MODEL pascalTriangle." No user should
292 instance it by itself. *)
293
294 col[set] IS_A factor;
295
296 END pTrow;
297
298 (* ---------------------------------------------------------- *)
299
300 MODEL pascalTriangle(
301 dim IS_A integer_constant;
302 ) REFINES ivpBase;
303
304 (* Creates a Pascal Triangle of dimension "dim."
305 A Pascal triangle gives the coefficients for the
306 polynomials resulting from expanding (x-y)^q,
307 for q = 0..dim. For dim=5, the triangle we form
308 here is
309
310 row(1): 1 1 1 1 1 1
311 row(2): 1 2 3 4 5
312 row(3): 1 3 6 10
313 row(4): 1 4 10
314 row(5): 1 5
315 row(6): 1
316
317 The coefficients for q=4 for (x-y)^4 are then in
318 the next to last place in each row, namely
319 1, 4, 6, 4, 1
320 as in
321 (1)x^4 - (4)x^3y + (6)x^2y^2 - (4)xy^3 + (1)y^4.
322 *)
323
324 FOR i IN [0..dim] CREATE
325 row[i] IS_A pTrow([0..dim-i]);
326 END FOR;
327
328 METHODS
329
330 (* ----- pascalTriangle ----- *)
331
332 METHOD default_self;
333 FOR i IN [0..dim] DO
334 row[0].col[i] := 1;
335 row[i].col[0] := 1;
336 END FOR;
337 FOR i IN [1..dim] DO
338 FOR j IN [1..dim-i] DO
339 row[i].col[j] := row[i-1].col[j]+row[i].col[j-1];
340 END FOR;
341 END FOR;
342 END default_self;
343
344 (* ----- pascalTriangle ----- *)
345
346 METHOD specify;
347 FOR i IN [0..dim] DO
348 FOR j IN [0..dim-i] DO
349 row[i].col[j].fixed := TRUE;
350 END FOR;
351 END FOR;
352 END specify;
353
354 END pascalTriangle;
355
356 (* ---------------------------------------------------------- *)
357
358 MODEL integrationPoint REFINES ivpBase;
359
360 (* ---------------------------------------------------------- *)
361 (* ------ the user models should be based on this one ------- *)
362 (* ---------------------------------------------------------- *)
363
364 (* this model provides storage for the predicted variables,
365 y, and the state derivatives, dydx, at a point related
366 to the independent variable x. NOTE: y, dydx, AND x MUST
367 BE DIMENSIONLESS VARIABLES.
368
369 nStates is the total number of states for the problem.
370 nPredVars is the total of state and algebraic variables
371 that the system should predict when stepping to the next
372 point. If there are no predicted algebraic variables,
373 set it equal to nStates. You must map the predicted
374 algebraic variables into y[nState+1..nPredVars].
375
376 dx is this x less the x for the current point.
377
378 The user should write two forms of a model that refines
379 this one to write the equations for the physical model for
380 the problem.
381
382 The first form should establish the types for the problem
383 variables by using IS_REFINED_TO statements for x, all y,
384 and all dydt variables. This form will be used for all
385 points except the current point.
386
387 The second form should refine the first form and add in the
388 model equations - both the dynamic and the algebraic
389 equations. This form will be for the current point in
390 the model. *)
391
392 (* ------------------------------------------------------
393 -------------- test for index problems ----------------
394 -------------------------------------------------------
395
396 The user can use (possibly by refining) the method
397 provided below called testForIndexProblem to test the
398 problem to see if it has an index greater than one. This
399 method should set the degrees of freedom for the currentPoint
400 to fix the independent variable x and the states y and to
401 calculate the time derivatives dydx and all the algebraic
402 variables. The user should then send only the currentPoint
403 instance to the solver. If ASCEND reports the instance for
404 the currentPoint to be structurally or numerically singular
405 (use the tool Solver/Analyze/Find dependent eqns), the
406 problem has an index greater than one and must be
407 reformulated as the error control of these numerical
408 methods will be inadequate. *)
409
410 nStates IS_A integer_constant;
411 nPredVars IS_A integer_constant;
412 x IS_A factor;
413 y[1..nPredVars] IS_A factor;
414 dydx[1..nStates] IS_A factor;
415 dx IS_A factor;
416
417 METHODS
418
419 (* ----- integrationPoint ----- *)
420
421 METHOD default_self;
422 x := 0.0;
423 y[1..nPredVars] := 0.0;
424 dydx[1..nStates] := 0.0;
425 dx := 0.0;
426 END default_self;
427
428 (* ----- integrationPoint ----- *)
429
430 METHOD specify;
431 x.fixed := TRUE;
432 y[1..nPredVars].fixed := TRUE;
433 dydx[1..nStates].fixed := TRUE;
434 dx.fixed := TRUE;
435 END specify;
436
437 (* ----- integrationPoint ----- *)
438
439 METHOD testForIndexProblem;
440 RUN specify;
441 dydx[1..nStates].fixed := FALSE;
442 y[nStates+1..nPredVars].fixed := FALSE;
443 END testForIndexProblem;
444
445 END integrationPoint;
446
447 (* ---------------------------------------------------------- *)
448
449 MODEL coefVect(
450 nPredVars WILL_BE integer_constant;
451 ) REFINES ivpBase;
452
453 (* a vector for storing the polynomial coefficient for all
454 state variables *)
455
456 var[1..nPredVars] IS_A factor;
457
458 METHODS
459
460 (* ----- coefVect ----- *)
461
462 METHOD default_self;
463 var[1..nPredVars] := 0.0;
464 END default_self;
465
466 (* ----- coefVect ----- *)
467
468 METHOD specify;
469 var[1..nPredVars].fixed := TRUE;
470 END specify;
471
472 END coefVect;
473
474 (* ---------------------------------------------------------- *)
475
476 MODEL polyCoefs(
477 nPredVars WILL_BE integer_constant;
478 ) REFINES ivpBase;
479
480 (* storage for the all polynomial coefficients for the numerical
481 method. Each predicted variable has one polynomial stored
482 for it. Storage is in the form: a[j].var[i] for
483 coefficient j of the polynomial for variable i. *)
484
485 a[0..6] IS_A coefVect(nPredVars);
486
487 METHODS
488
489 (* ----- polyCoefs ----- *)
490
491 METHOD default_self;
492 RUN a[0..6].default_self;
493 END default_self;
494
495 (* ----- polyCoefs ----- *)
496
497 METHOD specify;
498 RUN a[0..6].specify;
499 END specify;
500
501 END polyCoefs;
502
503 (* ---------------------------------------------------------- *)
504
505 MODEL poly1Value(
506 iP WILL_BE integrationPoint;
507 p WILL_BE polyCoefs;
508 set WILL_BE set OF integer_constant;
509 ) REFINES ivpBase;
510
511 (* a model that produces the equations for a first order
512 polynomial - written in the form of a Taylor Series *)
513
514 FOR i IN [set] CREATE
515 eqnValue[i]: iP.y[i] = p.a[0].var[i]
516 +iP.dx*(p.a[1].var[i]);
517 END FOR;
518
519 END poly1Value;
520
521 (* ---------------------------------------------------------- *)
522
523 MODEL poly1Deriv(
524 iP WILL_BE integrationPoint;
525 p WILL_BE polyCoefs;
526 set WILL_BE set OF integer_constant;
527 ) REFINES ivpBase;
528
529 (* a model that produces the equations for the derivative for a
530 first order derivative polynomial
531 - written in the form of a Taylor Series *)
532
533 FOR i IN [set] CREATE
534 eqnDeriv[i]: iP.dydx[i] = p.a[1].var[i];
535 END FOR;
536
537 END poly1Deriv;
538
539 (* ---------------------------------------------------------- *)
540
541 MODEL poly2Value(
542 iP WILL_BE integrationPoint;
543 p WILL_BE polyCoefs;
544 set WILL_BE set OF integer_constant;
545 ) REFINES ivpBase;
546
547 (* a model that produces the equations for a second order
548 polynomial - written in the form of a Taylor Series *)
549
550 FOR i IN [set] CREATE
551 eqnValue[i]: iP.y[i] = p.a[0].var[i]
552 +iP.dx*(p.a[1].var[i]
553 +iP.dx*(p.a[2].var[i]/2));
554 END FOR;
555
556 END poly2Value;
557
558 (* ---------------------------------------------------------- *)
559
560 MODEL poly2Deriv(
561 iP WILL_BE integrationPoint;
562 p WILL_BE polyCoefs;
563 set WILL_BE set OF integer_constant;
564 ) REFINES ivpBase;
565
566 (* a model that produces the equations for the derivative for a
567 second order derivative polynomial
568 - written in the form of a Taylor Series *)
569
570 FOR i IN [set] CREATE
571 eqnDeriv[i]: iP.dydx[i] = p.a[1].var[i]
572 +iP.dx*(p.a[2].var[i]);
573 END FOR;
574
575 END poly2Deriv;
576
577 (* ---------------------------------------------------------- *)
578
579 MODEL poly3Value(
580 iP WILL_BE integrationPoint;
581 p WILL_BE polyCoefs;
582 set WILL_BE set OF integer_constant;
583 ) REFINES ivpBase;
584
585 (* a model that produces the equations for a third order
586 polynomial - written in the form of a Taylor Series *)
587
588 FOR i IN [set] CREATE
589 eqnValue[i]: iP.y[i] = p.a[0].var[i]
590 +iP.dx*(p.a[1].var[i]
591 +iP.dx*(p.a[2].var[i]/2
592 +iP.dx*(p.a[3].var[i]/6)));
593 END FOR;
594
595 END poly3Value;
596
597 (* ---------------------------------------------------------- *)
598
599 MODEL poly3Deriv(
600 iP WILL_BE integrationPoint;
601 p WILL_BE polyCoefs;
602 set WILL_BE set OF integer_constant;
603 ) REFINES ivpBase;
604
605 (* a model that produces the equations for the derivative for a
606 third order derivative polynomial
607 - written in the form of a Taylor Series *)
608
609 FOR i IN [set] CREATE
610 eqnDeriv[i]: iP.dydx[i] = p.a[1].var[i]
611 +iP.dx*(p.a[2].var[i]
612 +iP.dx*(p.a[3].var[i]/2));
613 END FOR;
614
615 END poly3Deriv;
616
617 (* ---------------------------------------------------------- *)
618
619 MODEL poly4Value(
620 iP WILL_BE integrationPoint;
621 p WILL_BE polyCoefs;
622 set WILL_BE set OF integer_constant;
623 ) REFINES ivpBase;
624
625 (* a model that produces the equations for a fourth order
626 polynomial - written in the form of a Taylor Series *)
627
628 FOR i IN [set] CREATE
629 eqnValue[i]: iP.y[i] = p.a[0].var[i]
630 +iP.dx*(p.a[1].var[i]
631 +iP.dx*(p.a[2].var[i]/2
632 +iP.dx*(p.a[3].var[i]/6
633 +iP.dx*(p.a[4].var[i]/24))));
634 END FOR;
635
636 END poly4Value;
637
638 (* ---------------------------------------------------------- *)
639
640 MODEL poly4Deriv(
641 iP WILL_BE integrationPoint;
642 p WILL_BE polyCoefs;
643 set WILL_BE set OF integer_constant;
644 ) REFINES ivpBase;
645
646 (* a model that produces the equations for the derivative for a
647 fourth order derivative polynomial
648 - written in the form of a Taylor Series *)
649
650 FOR i IN [set] CREATE
651 eqnDeriv[i]: iP.dydx[i] = p.a[1].var[i]
652 +iP.dx*(p.a[2].var[i]
653 +iP.dx*(p.a[3].var[i]/2
654 +iP.dx*(p.a[4].var[i]/6)));
655 END FOR;
656
657 END poly4Deriv;
658
659 (* ---------------------------------------------------------- *)
660
661 MODEL poly5Value(
662 iP WILL_BE integrationPoint;
663 p WILL_BE polyCoefs;
664 set WILL_BE set OF integer_constant;
665 ) REFINES ivpBase;
666
667 (* a model that produces the equations for a fifth order
668 polynomial - written in the form of a Taylor Series *)
669
670 FOR i IN [set] CREATE
671 eqnValue[i]: iP.y[i] = p.a[0].var[i]
672 +iP.dx*(p.a[1].var[i]
673 +iP.dx*(p.a[2].var[i]/2
674 +iP.dx*(p.a[3].var[i]/6
675 +iP.dx*(p.a[4].var[i]/24
676 +iP.dx*(p.a[5].var[i]/120)))));
677 END FOR;
678
679 END poly5Value;
680
681 (* ---------------------------------------------------------- *)
682
683 MODEL poly5Deriv(
684 iP WILL_BE integrationPoint;
685 p WILL_BE polyCoefs;
686 set WILL_BE set OF integer_constant;
687 ) REFINES ivpBase;
688
689 (* a model that produces the equations for the derivative for a
690 fifth order derivative polynomial
691 - written in the form of a Taylor Series *)
692
693 FOR i IN [set] CREATE
694 eqnDeriv[i]: iP.dydx[i] = p.a[1].var[i]
695 +iP.dx*(p.a[2].var[i]
696 +iP.dx*(p.a[3].var[i]/2
697 +iP.dx*(p.a[4].var[i]/6
698 +iP.dx*(p.a[5].var[i]/24))));
699 END FOR;
700
701 END poly5Deriv;
702
703 (* ---------------------------------------------------------- *)
704
705 MODEL poly6Value(
706 iP WILL_BE integrationPoint;
707 p WILL_BE polyCoefs;
708 set WILL_BE set OF integer_constant;
709 ) REFINES ivpBase;
710
711 (* a model that produces the equations for a sixth order
712 polynomial - written in the form of a Taylor Series *)
713
714 FOR i IN [set] CREATE
715 eqnValue[i]: iP.y[i] = p.a[0].var[i]
716 +iP.dx*(p.a[1].var[i]
717 +iP.dx*(p.a[2].var[i]/2
718 +iP.dx*(p.a[3].var[i]/6
719 +iP.dx*(p.a[4].var[i]/24
720 +iP.dx*(p.a[5].var[i]/120
721 +iP.dx*(p.a[6].var[i]/720))))));
722 END FOR;
723
724 END poly6Value;
725
726 (* ---------------------------------------------------------- *)
727
728 MODEL poly6Deriv(
729 iP WILL_BE integrationPoint;
730 p WILL_BE polyCoefs;
731 set WILL_BE set OF integer_constant;
732 ) REFINES ivpBase;
733
734 (* a model that produces the equations for the derivative for a
735 sixth order derivative polynomial
736 - written in the form of a Taylor Series *)
737
738 FOR i IN [set] CREATE
739 eqnDeriv[i]: iP.dydx[i] = p.a[1].var[i]
740 +iP.dx*(p.a[2].var[i]
741 +iP.dx*(p.a[3].var[i]/2
742 +iP.dx*(p.a[4].var[i]/6
743 +iP.dx*(p.a[5].var[i]/24
744 +iP.dx*(p.a[6].var[i]/120)))));
745 END FOR;
746
747 END poly6Deriv;
748
749 (* ---------------------------------------------------------- *)
750
751 MODEL multistepEqns(
752 iP[0..7] WILL_BE integrationPoint;
753 p WILL_BE polyCoefs;
754 ) REFINES ivpBase;
755
756 (* ----- multistepEqns -----
757 Set aside space for a 6th order multistep integration method.
758
759 Complicating these models is the need for it to exist for
760 all orders q possible for the polynomials one is fitting.
761 For this model, q is to be selected interactively by the
762 person (or script) solving.
763
764 Define
765 (1) y(x) = poly(coef ai, indep var x, order q)
766 (2) f(y(x)) = derivative of poly(...) wrt x
767 where y(x) and f(y(x)) are values computed by the
768 model. *)
769
770 (* ----- multistepEqns -----
771 define the sets we will need *)
772
773 states,
774 nonStates IS_A set OF integer_constant;
775 states :== [1..iP[0].nStates];
776 nonStates :== [iP[0].nStates+1..iP[0].nPredVars];
777
778 (* ----- multistepEqns -----
779 set up polynomials for all points. When setting up
780 a model involving polynomials of order j, we need to
781 create equations of type pjDerivS for the current
782 point and j-1 past points. We also need equations of
783 the type pjValueS and pjValueNS for the current point
784 and j past points. *)
785
786
787 FOR q IN [0..0] CREATE
788 p1DerivS[q] IS_A poly1Deriv(iP[q], p, states);
789 END FOR;
790 FOR q IN [0..1] CREATE
791 p1ValueS[q] IS_A poly1Value(iP[q], p, states);
792 p1ValueNS[q] IS_A poly1Value(iP[q], p, nonStates);
793
794 p2DerivS[q] IS_A poly2Deriv(iP[q], p, states);
795 END FOR;
796 FOR q IN [0..2] CREATE
797 p2ValueS[q] IS_A poly2Value(iP[q], p, states);
798 p2ValueNS[q] IS_A poly2Value(iP[q], p, nonStates);
799
800 p3DerivS[q] IS_A poly3Deriv(iP[q], p, states);
801 END FOR;
802 FOR q IN [0..3] CREATE
803 p3ValueS[q] IS_A poly3Value(iP[q], p, states);
804 p3ValueNS[q] IS_A poly3Value(iP[q], p, nonStates);
805
806 p4DerivS[q] IS_A poly4Deriv(iP[q], p, states);
807 END FOR;
808 FOR q IN [0..4] CREATE
809 p4ValueS[q] IS_A poly4Value(iP[q], p, states);
810 p4ValueNS[q] IS_A poly4Value(iP[q], p, nonStates);
811
812 p5DerivS[q] IS_A poly5Deriv(iP[q], p, states);
813 END FOR;
814 FOR q IN [0..5] CREATE
815 p5ValueS[q] IS_A poly5Value(iP[q], p, states);
816 p5ValueNS[q] IS_A poly5Value(iP[q], p, nonStates);
817
818 p6DerivS[q] IS_A poly6Deriv(iP[q], p, states);
819 END FOR;
820 FOR q IN [0..6] CREATE
821 p6ValueNS[q] IS_A poly6Value(iP[q], p, nonStates);
822 p6ValueS[q] IS_A poly6Value(iP[q], p, states);
823 END FOR;
824
825 METHODS
826
827 (* ----- multistepEqns ----- *)
828
829 METHOD default_self;
830 END default_self;
831
832 (* ----- multistepEqns ----- *)
833
834 METHOD specify;
835 END specify;
836
837 END multistepEqns;
838
839 (* ---------------------------------------------------------- *)
840
841 MODEL ivpStep REFINES ivpBase;
842
843 (* ----- ivpStep -----
844 The user should implement a model that refines this model
845 and use it to provide the following methods
846
847 valuesInitialPt
848 specifyInitialPt
849 valuesStepping
850 specifyStepping
851
852 suitable for it. *)
853
854 (* ----- ivpStep -----
855 This model sets up the ivpStep model, allocating all the
856 space needed for all the integration points and polynomial
857 coefficients that the multistep methods need. It also
858 provides the following methods:
859
860 stepX
861 start
862 incrementPolyOrder
863 decrementPolyOrder
864 setMethodToBdf
865 setMethodToMs
866
867 stepX: prepares the model to be used to take the next
868 step in the independent variable. It shifts all the
869 state variables and their derivatives back one step.
870 The Taylor Series polynomial is always written around
871 the current point, which means we have to alter them
872 to account for their "zero" being shifted by deltaX,
873 the step about to be taken. This shifting is done
874 here. Finally, this method uses these shifted
875 polynomials to estimate the values for the states and
876 their derivatives at the new current time point.
877
878 start: should be run just following setting the values
879 for the initial point to start the model at the first
880 integration step it is to take.
881
882 incrementPolyOrder: increase the order for the
883 polynomial to be used for integrating, unless it is
884 already at the maximum order allowed - here 6.
885
886 decrementPolyOrder: decrease the order for the
887 polynomial to be used for integrating, unless the
888 order is already one, the least allowed.
889
890 setMethodToBdf: set the method to use to be bdf.
891
892 setMethodToAm: set the method to use to be MS.
893 *)
894
895 (* ----- ivpStep ----- *)
896
897 iP[0..7] IS_A integrationPoint;
898 nStates ALIASES iP[0].nStates;
899 nPredVars ALIASES iP[0].nPredVars;
900 p IS_A polyCoefs(nPredVars);
901 ms IS_A multistepEqns(iP, p);
902
903 usePolyOrder IS_A integer;
904 useMethod IS_A symbol;
905 maxStep, maxStepNew IS_A factor;
906 maxStepOrder IS_A symbol;
907
908
909 (* ----- ivpStep -----
910 Define
911 (1) y(x) = poly(coef ai, indep var x, order q)
912 (2) f(y(x)) = derivative of poly(...) wrt x
913 where y(x) and f(y(x)) are values computed by the
914 model.
915
916 *)
917
918 (* ----- ivpStep -----
919 Create the bdf model for stiff problems.
920
921 CASE 'bdf',j: is for writing a bdf method based
922 on j-th order polynomials.
923
924 Write (1) for all points. In addition write
925 (2) at the current point for the states.
926 *)
927
928 WHEN (useMethod, usePolyOrder)
929
930 CASE 'bdf',1:
931 FOR q IN [0..1] CREATE
932 USE ms.p1ValueS[q];
933 USE ms.p1ValueNS[q];
934 END FOR;
935 USE ms.p1DerivS[0];
936 CASE 'bdf',2:
937 FOR q IN [0..2] CREATE
938 USE ms.p2ValueS[q];
939 USE ms.p2ValueNS[q];
940 END FOR;
941 USE ms.p2DerivS[0];
942 CASE 'bdf',3:
943 FOR q IN [0..3] CREATE
944 USE ms.p3ValueS[q];
945 USE ms.p3ValueNS[q];
946 END FOR;
947 USE ms.p3DerivS[0];
948 CASE 'bdf',4:
949 FOR q IN [0..4] CREATE
950 USE ms.p4ValueS[q];
951 USE ms.p4ValueNS[q];
952 END FOR;
953 USE ms.p4DerivS[0];
954 CASE 'bdf',5:
955 FOR q IN [0..5] CREATE
956 USE ms.p5ValueS[q];
957 USE ms.p5ValueNS[q];
958 END FOR;
959 USE ms.p5DerivS[0];
960 CASE 'bdf',6:
961 FOR q IN [0..6] CREATE
962 USE ms.p6ValueS[q];
963 USE ms.p6ValueNS[q];
964 END FOR;
965 USE ms.p6DerivS[0];
966
967 (* ----- ivpStep -----
968 Create the Adams Moulton model for non-stiff problem.
969
970 For the states using the Adams Moulton method, write
971 (1) for the current [0] and previous point [1]. Write
972 (2) for all points [0..order-1].
973
974 For predicted non-state variables, write (1) for
975 all points.
976 *)
977
978 CASE 'am',1:
979 FOR q IN [0..1] CREATE
980 USE ms.p1ValueS[q];
981 END FOR;
982 FOR q IN [0..0] CREATE
983 USE ms.p1DerivS[q];
984 USE ms.p1ValueNS[q];
985 END FOR;
986 USE ms.p1ValueNS[1];
987 CASE 'am',2:
988 FOR q IN [0..1] CREATE
989 USE ms.p2ValueS[q];
990 END FOR;
991 FOR q IN [0..1] CREATE
992 USE ms.p2DerivS[q];
993 USE ms.p2ValueNS[q];
994 END FOR;
995 USE ms.p2ValueNS[2];
996 CASE 'am',3:
997 FOR q IN [0..1] CREATE
998 USE ms.p3ValueS[q];
999 END FOR;
1000 FOR q IN [0..2] CREATE
1001 USE ms.p3DerivS[q];
1002 USE ms.p3ValueNS[q];
1003 END FOR;
1004 USE ms.p3ValueNS[3];
1005 CASE 'am',4:
1006 FOR q IN [0..1] CREATE
1007 USE ms.p4ValueS[q];
1008 END FOR;
1009 FOR q IN [0..3] CREATE
1010 USE ms.p4DerivS[q];
1011 USE ms.p4ValueNS[q];
1012 END FOR;
1013 USE ms.p4ValueNS[4];
1014 CASE 'am',5:
1015 FOR q IN [0..1] CREATE
1016 USE ms.p5ValueS[q];
1017 END FOR;
1018 FOR q IN [0..4] CREATE
1019 USE ms.p5DerivS[q];
1020 USE ms.p5ValueNS[q];
1021 END FOR;
1022 USE ms.p5ValueNS[5];
1023 CASE 'am',6:
1024 FOR q IN [0..1] CREATE
1025 USE ms.p6ValueS[q];
1026 END FOR;
1027 FOR q IN [0..5] CREATE
1028 USE ms.p6DerivS[q];
1029 USE ms.p6ValueNS[q];
1030 END FOR;
1031 USE ms.p6ValueNS[6];
1032 END WHEN;
1033
1034 currentPt ALIASES iP[0];
1035 f IS_A factorial(6);
1036 pT IS_A pascalTriangle(6);
1037 deltaX IS_A factor;
1038 maxDeltaX IS_A factor;
1039 stopX IS_A factor;
1040 dxPow[0..6] IS_A factor;
1041 dxMax[['down','same','up']][1..nPredVars] IS_A factor;
1042 maxNominalSteppingError IS_A factor;
1043 lastHighestDeriv[1..nPredVars] IS_A factor;
1044
1045 stopConds IS_A set OF integer_constant;
1046
1047 newStopCondHit,
1048 stopCondHit,
1049 stopXHit,
1050 thisIsTheFinalStep IS_A integer;
1051
1052 eps IS_A factor;
1053
1054 stepCalcEqn: iP[0].x = iP[1].x + deltaX + iP[0].dx;
1055
1056 METHODS
1057
1058 (* ----- ivpStep ----- *)
1059
1060 METHOD default_self;
1061 RUN iP[0..7].default_self;
1062 RUN p.default_self;
1063 RUN ms.default_self;
1064 RUN f.default_self;
1065 RUN pT.default_self;
1066 eps := 1.0e-6;
1067 END default_self;
1068
1069 (* ----- ivpStep ----- *)
1070
1071 METHOD specify;
1072 RUN iP[0..7].specify;
1073 RUN p.specify;
1074 RUN ms.specify;
1075 RUN f.specify;
1076 RUN pT.specify;
1077 deltaX.fixed := TRUE;
1078 dxPow[0..6].fixed := TRUE;
1079 iP[0].y[1..nPredVars].fixed := FALSE;
1080 iP[0].dydx[1..nStates].fixed := FALSE;
1081 p.a[0..usePolyOrder].var[1..nPredVars].fixed := FALSE;
1082 iP[0..usePolyOrder].dx.fixed := TRUE;
1083 iP[0].x.fixed := FALSE;
1084 dxMax[['down','same','up']][1..nPredVars].fixed := TRUE;
1085 maxNominalSteppingError.fixed := TRUE;
1086 lastHighestDeriv[1..nPredVars].fixed := TRUE;
1087 eps.fixed := TRUE;
1088 END specify;
1089
1090 (* ----- ivpStep ----- *)
1091
1092 METHOD values;
1093 (* ----- ivpStep -----
1094 run this method after establishing the initial
1095 conditions for the model *)
1096
1097 newStopCondHit := 0;
1098 stopCondHit := 0;
1099 stopXHit := 0;
1100 thisIsTheFinalStep := 0;
1101
1102 FOR j IN [1..6] DO
1103 iP[j].x := 0;
1104 iP[j].y[1..nPredVars] := 0;
1105 iP[j].dydx[1..nStates] := 0;
1106 END FOR;
1107
1108 usePolyOrder := 1;
1109 FOR i IN [1..nPredVars] DO
1110 p.a[0].var[i] := iP[0].y[i];
1111 END FOR;
1112 FOR i IN [1..nStates] DO
1113 p.a[1].var[i] := iP[0].dydx[i];
1114 END FOR;
1115 END values;
1116
1117 (* ----- ivpStep ----- *)
1118
1119 METHOD incrementUsePolyOrder;
1120 IF usePolyOrder <= 5 THEN
1121 usePolyOrder := usePolyOrder +1;
1122 p.a[usePolyOrder].var[1..nPredVars].fixed := FALSE;
1123 END IF;
1124 END incrementUsePolyOrder;
1125
1126 (* ----- ivpStep ----- *)
1127
1128 METHOD decrementUsePolyOrder;
1129 IF usePolyOrder >= 2 THEN
1130 p.a[usePolyOrder].var[1..nPredVars] := 0.0;
1131 p.a[usePolyOrder].var[1..nPredVars].fixed := TRUE;
1132 usePolyOrder := usePolyOrder -1;
1133 END IF;
1134 END decrementUsePolyOrder;
1135
1136 (* ----- ivpStep ----- *)
1137
1138 METHOD setUseMethodToBdf;
1139 useMethod := 'bdf';
1140 END setUseMethodToBdf;
1141
1142 (* ----- ivpStep ----- *)
1143
1144 METHOD setUseMethodToAm;
1145 useMethod := 'am';
1146 END setUseMethodToAm;
1147
1148 (* ----- ivpStep ----- *)
1149
1150 METHOD stepX;
1151 (* ----- ivpStep -----
1152 run stepX to increment the independent variable and estimate
1153 the values for the states and predicted algebraic variables
1154 at that point *)
1155
1156 (* ----- ivpStep -----
1157 remember the last highest derivative, as the difference in
1158 it this time and last allows us to estimate the derivative
1159 that is one higher than we are actually computing. We need
1160 this derivative estimate for step size control. We use
1161 nominal values to make this estimate dimensionless. *)
1162
1163 FOR i IN [1..nPredVars] DO
1164 lastHighestDeriv[i] := p.a[usePolyOrder].var[i];
1165 END FOR;
1166
1167 (* ----- ivpStep -----
1168 we write the polynomials for the multistep method(s) as
1169 Taylor Series expansions about the current point, which is
1170 the point in iP[0]. When we step to a new current point,
1171 we must shift the polynomial coefficients accordingly so the
1172 shifted polynomial predicts the shifted past points correctly.
1173 Note, that after we shift the polnomial to be a Taylor series
1174 expansion around the new current point, the leading
1175 coefficient (a[0]) for each polynomial is the value
1176 for the corresponding state variable at the current point so
1177 this shifting is in fact a prediction method for the
1178 multistep process.
1179 *)
1180
1181 (* ----- ivpStep -----
1182 shift the predicted variables and dydx information and the
1183 time *)
1184
1185 FOR j IN [0..usePolyOrder] DO
1186 iP[usePolyOrder-j+1].x := iP[usePolyOrder-j].x;
1187 FOR i IN [1..nPredVars] DO
1188 iP[usePolyOrder-j+1].y[i] := iP[usePolyOrder-j].y[i];
1189 END FOR; (* i *)
1190 FOR i IN [1..nStates] DO
1191 iP[usePolyOrder-j+1].dydx[i] := iP[usePolyOrder-j].dydx[i];
1192 END FOR; (* i *)
1193 END FOR; (* j *)
1194
1195 iP[0].x := iP[1].x + deltaX;
1196 IF iP[0].x > stopX THEN
1197 iP[0].x := stopX;
1198 deltaX := iP[0].x - iP[1].x;
1199 END IF;
1200
1201 FOR j IN [0..usePolyOrder] DO
1202 iP[j].dx := iP[j].x - iP[0].x;
1203 END FOR; (* j *)
1204
1205 (* ----- ivpStep -----
1206 Estimate the polynomial coefficients for x shifted
1207 so new point is at x=0.
1208
1209 Update is based on Pascal's triangle and the
1210 shift, deltaX. The order we do this here allows
1211 us to do the update within the original data space
1212 for the polynomial.
1213 *)
1214
1215 FOR k IN [0..usePolyOrder] DO
1216 dxPow[0] := deltaX;
1217 FOR j IN [1..usePolyOrder-k] DO
1218 FOR i IN [1..nPredVars] DO
1219 p.a[k].var[i] := p.a[k].var[i]
1220 +dxPow[j-1]*pT.row[j].col[k]
1221 *p.a[k+j].var[i]*f.fac[k]/f.fac[j];
1222 END FOR;
1223 dxPow[j] := deltaX*dxPow[j-1];
1224 END FOR;
1225 END FOR;
1226
1227 (* ----- ivpStep -----
1228 The new coefficients a[0] provide the values for
1229 the new point and a[1] their derivatives. *)
1230
1231 FOR i IN [1..nPredVars] DO
1232 iP[0].y[i] := p.a[0].var[i];
1233 END FOR;
1234 FOR i IN [1..nStates] DO
1235 iP[0].dydx[i] := p.a[1].var[i];
1236 END FOR;
1237
1238 END stepX;
1239
1240
1241 METHOD computeMaxNominalStepsForEachVariable;
1242
1243 (* ----- ivpStep -----
1244 compute the maximum step size one can use for each state and
1245 predicted algebraic variable. The highest order term in each
1246 Taylor series expansion estimates error.
1247
1248 Assume q is the current model order. We shall do this
1249 prediction for three cases: if we were to reduce the polynomial
1250 order to q-1, leave it at q, and increase it to q+1. (Yes, the
1251 step size may be larger for a lower polynomial order.)
1252
1253 We use nominal values to keep a correct dimensionality in all
1254 equations. *)
1255
1256 (* ----- ivpStep -----
1257 when reducing the polynomial order to q-1, we assume the Taylor
1258 series coefficient one less than we are now using stays constant
1259 *)
1260
1261 maxStepOrder := 'down';
1262
1263 maxStep := stopX - iP[1].x;
1264 IF maxStep >= maxDeltaX THEN
1265 maxStep := maxDeltaX;
1266 END IF;
1267
1268 IF usePolyOrder >=2 THEN
1269 FOR i IN [1..nPredVars] DO
1270 (*
1271 dxMax['down'][i] := currentPt.x.nominal
1272 *(abs(f.fac[usePolyOrder-1]
1273 *(maxNominalSteppingError*currentPt.y[i].nominal
1274 /currentPt.x.nominal^(usePolyOrder-1))
1275 /p.a[usePolyOrder-1].var[i]))
1276 ^(1.0/(usePolyOrder-1));
1277 *)
1278
1279 dxMax['down'][i] := (abs(f.fac[usePolyOrder-1]
1280 *maxNominalSteppingError/p.a[usePolyOrder-1].var[i]))
1281 ^(1.0/(usePolyOrder-1));
1282
1283 IF dxMax['down'][i] < maxStep THEN
1284 maxStep := dxMax['down'][i];
1285 END IF;
1286 END FOR;
1287
1288 ELSE
1289 dxMax['down'][1..nPredVars] := 0.0;
1290 maxStep := 0.0;
1291 END IF;
1292
1293
1294 FOR i IN [2..nPredVars] DO
1295 END FOR;
1296
1297 (* ----- ivpStep ----- *)
1298
1299 maxStepNew := stopX - iP[1].x;
1300 IF maxStepNew >= maxDeltaX THEN
1301 maxStepNew := maxDeltaX;
1302 END IF;
1303
1304 FOR i IN [1..nPredVars] DO
1305
1306 dxMax['same'][i] := (abs(f.fac[usePolyOrder]
1307 *maxNominalSteppingError
1308 /p.a[usePolyOrder].var[i]))^(1.0/usePolyOrder);
1309
1310 (* dxMax['same'][i] := currentPt.x.nominal*(abs(f.fac[usePolyOrder]
1311 *(maxNominalSteppingError*currentPt.y[i].nominal
1312 /currentPt.x.nominal^usePolyOrder)
1313 /p.a[usePolyOrder].var[i]))^(1.0/usePolyOrder);*)
1314
1315 IF dxMax['same'][i] < maxStepNew THEN
1316 maxStepNew := dxMax['same'][i];
1317 END IF;
1318 END FOR;
1319
1320 (* Find if 'same' step is better for choosing next step size*)
1321
1322 IF maxStepNew > maxStep THEN
1323 maxStep := maxStepNew;
1324 maxStepOrder := 'same';
1325 END IF;
1326
1327 (* ----- ivpStep -----
1328 we estimate the next higher derivative as the difference between
1329 the the q-th derivatives at this and the previous point, divided
1330 by the step size (i.e., delta/delta x approximates d/dx) *)
1331
1332 IF usePolyOrder <= 5 THEN
1333 maxStepNew := stopX - iP[1].x;
1334 IF maxStepNew >= maxDeltaX THEN
1335 maxStepNew := maxDeltaX;
1336 END IF;
1337
1338 FOR i IN [1..nPredVars] DO
1339
1340 dxMax['up'][i] := (abs(f.fac[usePolyOrder+1]
1341 *maxNominalSteppingError
1342 *deltaX
1343 /(p.a[usePolyOrder].var[i]
1344 -lastHighestDeriv[i])))
1345 ^(1.0/(usePolyOrder+1));
1346
1347 (* dxMax['up'][i] := currentPt.x.nominal
1348 *(abs(f.fac[usePolyOrder+1]
1349 *(maxNominalSteppingError*currentPt.y[i].nominal
1350 /currentPt.x.nominal^(usePolyOrder+1))
1351 *deltaX
1352 /(p.a[usePolyOrder].var[i]
1353 -lastHighestDeriv[i]*currentPt.y[i].nominal
1354 /deltaX.nominal^usePolyOrder)))
1355 ^(1.0/(usePolyOrder+1)); *)
1356
1357 IF dxMax['up'][i] < maxStepNew THEN
1358 maxStepNew := dxMax['up'][i];
1359 END IF;
1360 END FOR;
1361
1362 (* Find if 'up' step is better for choosing next step size*)
1363
1364 IF maxStepNew > maxStep THEN
1365 maxStep := maxStepNew;
1366 maxStepOrder := 'up';
1367 END IF;
1368 ELSE
1369 dxMax['up'][1..nPredVars] := 0.0;
1370 END IF;
1371
1372 (* set step size and poly order for next step *)
1373
1374 deltaX := maxStep;
1375
1376 IF maxStepOrder = 'down' THEN
1377 RUN decrementUsePolyOrder;
1378 ELSE
1379 IF maxStepOrder = 'up' THEN
1380 RUN incrementUsePolyOrder;
1381 END IF;
1382 END IF;
1383
1384 END computeMaxNominalStepsForEachVariable;
1385
1386 METHOD setStopConditions;
1387
1388 newStopCondHit := 0;
1389 iP[0].x.lower_bound := iP[1].x.lower_bound;
1390 iP[0].x.upper_bound := iP[1].x.upper_bound;
1391
1392 FOR i IN stopConds DO
1393 iP[0].y[i].fixed := FALSE;
1394 IF (newStopCondHit == 0) AND (abs(iP[0].y[i]) >= eps*iP[0].y[i].nominal) THEN
1395 IF (iP[0].y[i]*iP[1].y[i] <= -(eps*iP[0].y[i].nominal)^2) THEN
1396 iP[0].y[i] := 0.0;
1397 iP[0].dx.lower_bound := iP[1].dx;
1398 iP[0].dx.upper_bound := 0.0;
1399 iP[0].y[i].fixed := TRUE;
1400 iP[0].dx.fixed := FALSE;
1401 newStopCondHit := 1;
1402 stopCondHit := i;
1403 END IF;
1404 END IF;
1405 END FOR;
1406
1407 IF newStopCondHit == 0 THEN
1408 IF maxStep >= (stopX - iP[1].x)*(1.0 - eps) THEN
1409 stopXHit := 1;
1410 END IF;
1411 IF (stopCondHit >= 1) OR (stopXHit == 1) THEN
1412 deltaX := iP[0].x - iP[1].x;
1413 thisIsTheFinalStep := 1;
1414 iP[0].dx.fixed := TRUE;
1415 END IF;
1416 END IF;
1417
1418 END setStopConditions;
1419
1420 END ivpStep;

john.pye@anu.edu.au
ViewVC Help
Powered by ViewVC 1.1.22