| 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; |
Parent Directory
| 
Revision Log