/[ascend]/trunk/models/bvp.a4l
ViewVC logotype

Contents of /trunk/models/bvp.a4l

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2651 - (show annotations) (download) (as text)
Thu Dec 13 07:29:48 2012 UTC (11 years, 6 months ago) by jpye
File MIME type: text/x-ascend
File size: 16958 byte(s)
Fixing GPL header, removing postal address (rpmlint incorrect-fsf-address)
1 (* ASCEND modelling environment
2 Copyright (C) 1997 Benjamin A Allan
3 Copyright (C) 1998, 2007 Carnegie Mellon University
4
5 This program is free software; you can redistribute it and/or modify
6 it under the terms of the GNU General Public License as published by
7 the Free Software Foundation; either version 2, or (at your option)
8 any later version.
9
10 This program is distributed in the hope that it will be useful,
11 but WITHOUT ANY WARRANTY; without even the implied warranty of
12 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 GNU General Public License for more details.
14
15 You should have received a copy of the GNU General Public License
16 along with this program. If not, see <http://www.gnu.org/licenses/>.
17 *)
18 REQUIRE "atoms.a4l";
19 (*
20 This file contains models that implement an integration scheme completely
21 within ASCEND as a system of equations. Because data points are all
22 stored in memory, this approach uses a lot more memory (although with
23 relation sharing, sparse matrix storage, etc, it's not completely terrible)
24
25 @TODO more documentation here.
26
27 @NOTE there is a test model implemented at the bottom of the file. See
28 'test_dynamics' and 'bvp_test'. Note also that the model 'plotbvp.a4c'
29 includes this test model and provides graphical output from the integration.
30
31 @NOTE some old comments were removed here as they confused the issue by
32 referring to the LSODE integration engine.
33
34 @NOTE The LSODE integration engine is provided as part of ASCEND and
35 offers a more efficient way to solve integration models. Note that you will
36 need to set up your model quite differently if you want to use that approach.
37
38 Original code 1991.
39 Decoupled from IVP code, rewritten in ASCEND IV -- Ben Allan, Mar 1997
40
41 by Peter Piela, Boyd T. Safrit, Joseph J. Zaher, 1991
42 remodelled by Benjamin A. Allan
43 *)
44
45 MODEL bvp_base REFINES cmumodel;
46 NOTES 'purpose' SELF {
47 This MODEL anchors the boundary value problem library.
48 This library allows you to solve differential/algebraic
49 equation systems using standard collocation methods.
50 See example bvp_test.
51 }
52 END NOTES;
53 END bvp_base;
54
55 (*-------------------------------
56 bvp_point: the collocation user MODEL interface
57 *)
58 MODEL bvp_point REFINES bvp_base;
59 NOTES 'purpose' SELF {
60 Generic definition of the differential equation system,
61 Dy[i] = dydx[i], i=1..n_eq, D = d/dx.
62 Create a refinement of this type that includes a part which
63 is your DAE MODEL and merge your variables with x,y,dydx as
64 needed. See MODEL bvp_test for an example.
65 Your refinement should define all the standard methods
66 default_self, etc.
67 }
68 'usage' SELF {
69 Note that the integration will behave much better if the
70 typical values of your state variables are all about the
71 same size. Integrators can lie if you give them poorly
72 scaled problems. Dimensionless variables are generally best
73 to use if you can easily formulate your MODEL in those
74 terms.
75 }
76 END NOTES;
77
78 n_eq "number of implicit ODE's" IS_A integer_constant;
79
80 dydx[1..n_eq] "derivative variables",
81 y[1..n_eq] "state variables",
82 x "independent variable"
83 IS_A solver_var;
84
85 METHODS
86 METHOD default_self;
87 (* do not assign anything as these are merged with
88 * user's variables which are all assigned.
89 *)
90 END default_self;
91
92 METHOD specify;
93 (* The state variables are set
94 * FIXED, while their derivatives are FREE.
95 *)
96 FIX x;
97 FIX y[1..n_eq];
98 FREE dydx[1..n_eq];
99 END specify;
100 END bvp_point;
101
102 (*-------------------------------
103 BVP internal models
104 *)
105
106 MODEL propagation(
107 n_point WILL_BE integer_constant;
108 eval[0..n_point] WILL_BE bvp_point;
109 n_eq WILL_BE integer_constant;
110 )WHERE(
111 n_eq == eval[0].n_eq;
112 )REFINES bvp_base;
113
114 (* eval[0..n_point] ARE_ALIKE;
115 * Eval[i] must be deeply alike, so use of ARE_ALIKE here is heuristic at best.
116 * Formal type compatibility does not ensure a well defined BVP.
117 * No longer an issue, since eval is coming in from outside where it
118 * was constructed all at the same time.
119 *)
120
121 NOTES
122 'purpose' SELF {
123 This MODEL does nothing by itself, and should never be created.
124 Only the refinements of this MODEL should be created.
125 If you need a new collocation formula, create it by refining
126 propagation and then update the SELECT statement in
127 the MODEL integration.
128 }
129 'usage' n_order {
130 n_order did nothing.
131 It was mainly a carrier of information for the user.
132 Specifically, it gives the order of the error of the collocation.
133 in terms of h.
134 this information should just be kept in NOTES instead, like so.
135 }
136 END NOTES;
137
138
139 h IS_A solver_var;
140 (* propagate the dimensionality of x to h, though the
141 h, eval[0].x ARE_ALIKE;
142 * ASCEND environment would do the same thing for us automagically.
143 * We actually couldn't care less about the type of h, so
144 * long as it is dimensionally consistent with the x.
145 * or we could do it with an assignment:
146 h := eval[0].x;
147 *)
148 initial ALIASES eval[0];
149 final ALIASES eval[n_point];
150 final.x = initial.x + h;
151
152 METHODS
153 METHOD specify;
154 RUN eval[0..n_point].specify;
155 FIX h;
156 FREE final.x;
157 FREE eval[1..n_point].y[1..n_eq];
158 FREE eval[1..n_point].x;
159 END specify;
160 END propagation;
161
162
163 (*-------------------------------
164 Euler method
165 *)
166 MODEL euler(
167 n_point WILL_BE integer_constant;
168 eval[0..n_point] WILL_BE bvp_point;
169 n_eq WILL_BE integer_constant;
170 ) WHERE (
171 n_eq == eval[0].n_eq;
172 n_point == 1;
173 ) REFINES propagation();
174
175 NOTES
176 'order of error' SELF {
177 The error in the euler method is O(h).
178 }
179 END NOTES;
180
181 FOR i IN [1..n_eq] CREATE
182 eval[1].y[i] = eval[0].y[i] + h*eval[0].dydx[i];
183 END FOR;
184
185 METHODS
186 (* 'specify' inherited from propagation is correct *)
187 END euler;
188
189 (*-------------------------------
190 Trapezoidal method
191 *)
192 MODEL trapezoid(
193 n_point WILL_BE integer_constant;
194 eval[0..n_point] WILL_BE bvp_point;
195 n_eq WILL_BE integer_constant;
196 ) WHERE (
197 n_eq == eval[0].n_eq;
198 n_point == 1;
199 ) REFINES propagation();
200
201 NOTES
202 'order of error' SELF {
203 The error in the trapezoidal rule is O(h^2).
204 }
205 END NOTES;
206
207 FOR i IN [1..n_eq] CREATE
208 eval[1].y[i] = eval[0].y[i] + 0.5*h*
209 (eval[0].dydx[i] + eval[1].dydx[i]);
210 END FOR;
211
212 METHODS
213 (* 'specify' inherited from propagation is correct *)
214 END trapezoid;
215
216 (*-------------------------------
217 Midpoint method
218 *)
219 MODEL midpoint(
220 n_point WILL_BE integer_constant;
221 eval[0..n_point] WILL_BE bvp_point;
222 n_eq WILL_BE integer_constant;
223 ) WHERE (
224 n_eq == eval[0].n_eq;
225 n_point == 2;
226 ) REFINES propagation();
227
228 NOTES
229 'order of error' SELF {
230 The error in the midpoint rule is O(h^2).
231 }
232 END NOTES;
233
234 eval[1].x = initial.x + (1.0/2.0)*h;
235 FOR i IN [1..n_eq] CREATE
236 eval[1].y[i] = eval[0].y[i] + (1.0/2.0)*h*
237 eval[0].dydx[i];
238 eval[2].y[i] = eval[0].y[i] + h*
239 eval[1].dydx[i];
240 END FOR;
241
242 METHODS
243 (* 'specify' inherited from propagation is correct *)
244 END midpoint;
245
246 (*-------------------------------
247 Simpson's method
248 *)
249 MODEL simpsons(
250 n_point WILL_BE integer_constant;
251 eval[0..n_point] WILL_BE bvp_point;
252 n_eq WILL_BE integer_constant;
253 ) WHERE (
254 n_eq == eval[0].n_eq;
255 n_point == 2;
256 ) REFINES propagation();
257
258 NOTES
259 'order of error' SELF {
260 The error in simpson's rule is O(h^4).
261 }
262 END NOTES;
263
264 eval[1].x = initial.x + 0.5*h;
265 FOR i IN [1..n_eq] CREATE
266 eval[1].y[i] = eval[0].y[i] + 0.5*h*
267 ((5.0/12.0)*eval[0].dydx[i] +
268 (8.0/12.0)*eval[1].dydx[i] -
269 (1.0/12.0)*eval[2].dydx[i]);
270 eval[2].y[i] = eval[0].y[i] + h*
271 ((1.0/6.0)*eval[0].dydx[i] +
272 (4.0/6.0)*eval[1].dydx[i] +
273 (1.0/6.0)*eval[2].dydx[i]);
274 END FOR;
275
276 METHODS
277 (* 'specify' inherited from propagation is correct *)
278 END simpsons;
279
280 (*-------------------------------
281 4-point Runge-Kutta-Gill method
282 *)
283 MODEL runge_kutta_gill(
284 n_point WILL_BE integer_constant;
285 eval[0..n_point] WILL_BE bvp_point;
286 n_eq WILL_BE integer_constant;
287 ) WHERE (
288 n_eq == eval[0].n_eq;
289 n_point == 4;
290 ) REFINES propagation();
291
292 NOTES
293 'order of error' SELF {
294 The error in the 4 point Runge-Kutta-Gill method is O(h^4).
295 }
296 END NOTES;
297
298 eval[1].x = eval[2].x;
299 eval[3].x = eval[4].x;
300 eval[1].x = initial.x + 0.5*h;
301 FOR i IN [1..n_eq] CREATE
302 eval[1].y[i] = eval[0].y[i] + 0.5*h*
303 eval[0].dydx[i];
304 eval[2].y[i] = eval[0].y[i] + 0.5*h*
305 (0.41421356*eval[0].dydx[i] +
306 0.58578644*eval[1].dydx[i]);
307 eval[3].y[i] = eval[0].y[i] + h*
308 (-0.707106781*eval[1].dydx[i] +
309 1.707106781*eval[2].dydx[i]);
310 eval[4].y[i] = eval[0].y[i] + h*
311 (0.166666667*eval[0].dydx[i] +
312 0.097631073*eval[1].dydx[i] +
313 0.569035590*eval[2].dydx[i] +
314 0.166666667*eval[3].dydx[i]);
315 END FOR;
316
317 METHODS
318 (* 'specify' inherited from propagation is correct *)
319 END runge_kutta_gill;
320
321 (*-------------------------------
322 4-point Gauss method
323 *)
324 MODEL gauss(
325 n_point WILL_BE integer_constant;
326 eval[0..n_point] WILL_BE bvp_point;
327 n_eq WILL_BE integer_constant;
328 ) WHERE (
329 n_eq == eval[0].n_eq;
330 n_point == 4;
331 ) REFINES propagation();
332
333 NOTES
334 'order of error' SELF {
335 The error in 4 point gauss method is O(h^6).
336 }
337 END NOTES;
338
339 eval[1].x = initial.x + 0.11270167*h;
340 eval[2].x = initial.x + 0.5*h;
341 eval[3].x = initial.x + 0.88729833*h;
342 FOR i IN [1..n_eq] CREATE
343 eval[1].y[i] = initial.y[i] + 0.11270167*h*
344 (1.23236*eval[1].dydx[i] -
345 0.31922*eval[2].dydx[i] +
346 0.08686*eval[3].dydx[i]);
347 eval[2].y[i] = initial.y[i] + 0.5*h*
348 (0.60053*eval[1].dydx[i] +
349 0.44444*eval[2].dydx[i] -
350 0.04497*eval[3].dydx[i]);
351 eval[3].y[i] = initial.y[i] + 0.88729833*h*
352 (0.30203*eval[1].dydx[i] +
353 0.54144*eval[2].dydx[i] +
354 0.15653*eval[3].dydx[i]);
355 eval[4].y[i] = initial.y[i] + h*
356 (0.27778*eval[1].dydx[i] +
357 0.44444*eval[2].dydx[i] +
358 0.27778*eval[3].dydx[i]);
359 END FOR;
360
361 METHODS
362 (* 'specify' inherited from propagation is correct *)
363 END gauss;
364
365 (*-------------------------------
366 Lobatto method
367 *)
368 MODEL lobatto(
369 n_point WILL_BE integer_constant;
370 eval[0..n_point] WILL_BE bvp_point;
371 n_eq WILL_BE integer_constant;
372 ) WHERE (
373 n_eq == eval[0].n_eq;
374 n_point == 4;
375 ) REFINES propagation();
376
377 NOTES
378 'order of error' SELF {
379 The error in the 4 lobatto point method is O(h^8).
380 }
381 END NOTES;
382
383 eval[1].x = initial.x + 0.11270167*h;
384 eval[2].x = initial.x + 0.5*h;
385 eval[3].x = initial.x + 0.88729833*h;
386 FOR i IN [1..n_eq] CREATE
387 eval[1].y[i] = initial.y[i] + 0.11270167*h*
388 (0.428010*eval[0].dydx[i] +
389 0.602339*eval[1].dydx[i] -
390 0.044301*eval[2].dydx[i] +
391 0.029587*eval[3].dydx[i] -
392 0.015635*eval[4].dydx[i]);
393 eval[2].y[i] = initial.y[i] + 0.5*h*
394 (-0.06250*eval[0].dydx[i] +
395 0.68122*eval[1].dydx[i] +
396 0.44444*eval[2].dydx[i] -
397 0.12566*eval[3].dydx[i] +
398 0.06250*eval[4].dydx[i]);
399 eval[3].y[i] = initial.y[i] + 0.88729833*h*
400 (0.001986*eval[0].dydx[i] +
401 0.309303*eval[1].dydx[i] +
402 0.506524*eval[2].dydx[i] +
403 0.236552*eval[3].dydx[i] -
404 0.054365*eval[4].dydx[i]);
405 eval[4].y[i] = initial.y[i] + h*
406 (0.277778*eval[1].dydx[i] +
407 0.444444*eval[2].dydx[i] +
408 0.277778*eval[3].dydx[i]);
409 END FOR;
410
411 METHODS
412 (* 'specify' inherited from propagation is correct *)
413 END lobatto;
414
415 (*-------------------------------
416 INTEGRATION model
417 *)
418 MODEL integration(
419 nstep IS_A integer_constant;
420 npoint IS_A integer_constant;
421 method IS_A symbol_constant;
422 grid[0..nstep*npoint] WILL_BE bvp_point;
423 stepsize WILL_BE solver_var;
424 n_eq IS_A integer_constant;
425 ) WHERE (
426 ((method == 'Euler') AND (npoint == 1)) OR
427 ((method == 'trapezoidal') AND (npoint == 1)) OR
428 ((method == 'midpoint') AND (npoint == 2)) OR
429 ((method == 'Simpsons') AND (npoint == 2)) OR
430 ((method == 'Runge-Kutta-Gill') AND (npoint == 4)) OR
431 ((method == 'Gauss') AND (npoint == 4)) OR
432 ((method == 'Lobatto') AND (npoint == 4));
433 ) REFINES bvp_base;
434
435 NOTES 'usage' SELF {
436 This writes the equations for a mesh with nstep collocation
437 intervals. In each interval there npoint nodes.
438 The total mesh has (nstep * npoint + 1) nodes.
439
440 In order to use this MODEL, create a refinement of bvp_point
441 that contains the MODEL you want to integrate, then create
442 an array of them (grid) and pass it in along with the other options.
443 See the MODEL bvp_point for details of how to do this.
444 See the MODELs bvp_test and bvp_dynamics for an example.
445 }
446 END NOTES;
447 (* cook up the little arrays for propagation out of the
448 * big array from the user.
449 *)
450 FOR i IN [0..nstep-1] CREATE
451 step_nodes[i+1][pset[i]] ALIASES
452 (grid[(i*npoint)..((i+1)*npoint)])
453 WHERE pset[i] IS_A set OF integer_constant WITH_VALUE
454 (0..npoint);
455 END FOR;
456
457 (* merge all the sets. we need a better way of
458 * getting at the set defined for an array.
459 * we could get by without the merge.
460 *)
461 pset[0..nstep-1] ARE_THE_SAME;
462
463 (* create the generic step stuctures and refine to rule requested *)
464 FOR i IN [1..nstep] CREATE
465 step[i] IS_A propagation(npoint,step_nodes[i],n_eq);
466 END FOR;
467
468 (* refine to rule *)
469 SELECT (method)
470 CASE 'Euler':
471 FOR i IN [1..nstep] CREATE
472 step[i] IS_REFINED_TO
473 euler(npoint,step_nodes[i],n_eq);
474 END FOR;
475 CASE 'trapezoidal':
476 FOR i IN [1..nstep] CREATE
477 step[i] IS_REFINED_TO
478 trapezoid(npoint,step_nodes[i],n_eq);
479 END FOR;
480 CASE 'midpoint':
481 FOR i IN [1..nstep] CREATE
482 step[i] IS_REFINED_TO
483 midpoint(npoint,step_nodes[i],n_eq);
484 END FOR;
485 CASE 'Simpsons':
486 FOR i IN [1..nstep] CREATE
487 step[i] IS_REFINED_TO
488 simpsons(npoint,step_nodes[i],n_eq);
489 END FOR;
490 CASE 'Runge-Kutta-Gill':
491 FOR i IN [1..nstep] CREATE
492 step[i] IS_REFINED_TO
493 runge_kutta_gill(npoint,step_nodes[i],n_eq);
494 END FOR;
495 CASE 'Gauss':
496 FOR i IN [1..nstep] CREATE
497 step[i] IS_REFINED_TO
498 gauss(npoint,step_nodes[i],n_eq);
499 END FOR;
500 CASE 'Lobatto':
501 FOR i IN [1..nstep] CREATE
502 step[i] IS_REFINED_TO
503 lobatto(npoint,step_nodes[i],n_eq);
504 END FOR;
505 END SELECT;
506
507
508 METHODS
509 METHOD default_self;
510 RUN values;
511 END default_self;
512 METHOD values;
513 step[1..nstep].h := stepsize;
514 END values;
515 METHOD specify;
516 FIX stepsize;
517 FOR i IN [1..nstep] DECREASING DO
518 RUN step[i].specify;
519 END FOR;
520 END specify;
521 END integration;
522
523 (*-------------------------------
524 Demo and test models
525 *)
526
527 MODEL testbvp_base REFINES testcmumodel;
528 END testbvp_base;
529
530 MODEL test_dynamics REFINES testbvp_base;
531 x IS_A factor;
532 t IS_A time;
533 xdot IS_A frequency;
534 xdot = 2{1/s} * x * exp(-2{1/s}*t) * cos(t*1{rad/s});
535 METHODS
536 METHOD default_self;
537 RUN values;
538 END default_self;
539 METHOD specify;
540 FIX x;
541 FREE xdot;
542 FIX t;
543 END specify;
544 METHOD values;
545 xdot.lower_bound := -1e6{1/s};
546 END values;
547 END test_dynamics;
548
549 (*-------------------------------
550 'dynamic point'
551 *)
552 MODEL dynamic_point REFINES bvp_point;
553 n_eq :== 1; (* you can have as many as you like *)
554
555 (* replace the following line with an IS_A of your MODEL *)
556 your_model IS_A test_dynamics;
557
558 (* connect independent variables. If your MODEL is autonomous,
559 * you don't need to, but then you must call
560 * bvp_point::clear and bvp_point::specify in your methods
561 * to get x set properly.
562 *)
563 x, your_model.t ARE_THE_SAME;
564
565 (* connect your variables to the integration state variables *)
566 y[1], your_model.x ARE_THE_SAME;
567
568 (* connect your variables to the integration derivatives *)
569 dydx[1], your_model.xdot ARE_THE_SAME;
570
571 METHODS
572 METHOD specify;
573 RUN your_model.specify;
574 (* If you write your_model correctly, you don't need
575 to run bvp_point::specify because you've wired all
576 the bvp variables into some of yours.
577 (Not TRUE if your MODEL is autonomous, though).
578 RUN bvp_point::specify;
579 *)
580 END specify;
581 END dynamic_point;
582
583
584 MODEL bvp_test REFINES testbvp_base;
585 stepsize IS_A time;
586 npoint, nstep IS_A integer_constant;
587 npoint :== 1;
588 nstep :== 100;
589 nodes[0..nstep*npoint] IS_A dynamic_point;
590 initial ALIASES nodes[0];
591 final ALIASES nodes[nstep*npoint];
592 integral IS_A
593 integration(nstep,npoint,'Euler',nodes,stepsize,initial.n_eq);
594
595 (* we don't tie the test MODEL to the ASCEND plot MODEL *)
596
597 METHODS
598 METHOD default_self;
599 stepsize := 0.05{s};
600 initial.x := 0{s};
601 initial.y[1] := 1;
602 RUN nodes[0..nstep*npoint].your_model.default_self;
603 RUN integral.default_self;
604 END default_self;
605 METHOD specify;
606 FIX stepsize;
607 RUN integral.specify;
608 END specify;
609 METHOD values;
610 stepsize := 0.05{s};
611 initial.x := 0{s};
612 initial.y[1] := 1;
613 RUN nodes[0..nstep*npoint].your_model.values;
614 RUN integral.values;
615 END values;
616 METHOD on_load;
617 RUN reset; RUN default_self; RUN values;
618 END on_load;
619 METHOD self_test;
620 (* no tests yet *)
621 END self_test;
622 END bvp_test;
623 (* :ex: set ts=4: *)

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