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

Contents of /trunk/models/bvp.a4l

Parent Directory Parent Directory | Revision Log Revision Log


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

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