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

Annotation of /trunk/models/bvp.a4l

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2651 - (hide 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 johnpye 1158 (* 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 jpye 2651 along with this program. If not, see <http://www.gnu.org/licenses/>.
17 johnpye 1158 *)
18 aw0a 1 REQUIRE "atoms.a4l";
19     (*
20 johnpye 1158 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 aw0a 1
25 johnpye 1158 @TODO more documentation here.
26 aw0a 1
27 johnpye 1158 @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 aw0a 1
31 jpye 2314 @NOTE some old comments were removed here as they confused the issue by
32 johnpye 1158 referring to the LSODE integration engine.
33 aw0a 1
34 johnpye 1158 @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 aw0a 1
38 johnpye 1158 Original code 1991.
39     Decoupled from IVP code, rewritten in ASCEND IV -- Ben Allan, Mar 1997
40 aw0a 1
41 johnpye 1158 by Peter Piela, Boyd T. Safrit, Joseph J. Zaher, 1991
42     remodelled by Benjamin A. Allan
43 aw0a 1 *)
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 johnpye 1158
55     (*-------------------------------
56     bvp_point: the collocation user MODEL interface
57     *)
58 aw0a 1 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 johnpye 1158 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 aw0a 1
92 johnpye 1158 METHOD specify;
93     (* The state variables are set
94 jpye 2314 * FIXED, while their derivatives are FREE.
95 johnpye 1158 *)
96     FIX x;
97     FIX y[1..n_eq];
98     FREE dydx[1..n_eq];
99     END specify;
100 aw0a 1 END bvp_point;
101    
102 johnpye 1158 (*-------------------------------
103     BVP internal models
104     *)
105 aw0a 1
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 johnpye 1158 )WHERE(
111 aw0a 1 n_eq == eval[0].n_eq;
112 johnpye 1158 )REFINES bvp_base;
113    
114 aw0a 1 (* 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 johnpye 1158 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 aw0a 1 END propagation;
161    
162 johnpye 1158
163     (*-------------------------------
164     Euler method
165     *)
166 aw0a 1 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 johnpye 1158 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 aw0a 1
185     METHODS
186 johnpye 1158 (* 'specify' inherited from propagation is correct *)
187 aw0a 1 END euler;
188    
189 johnpye 1158 (*-------------------------------
190     Trapezoidal method
191     *)
192 aw0a 1 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 johnpye 1158 NOTES
202     'order of error' SELF {
203     The error in the trapezoidal rule is O(h^2).
204     }
205     END NOTES;
206 aw0a 1
207 johnpye 1158 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 aw0a 1
212 johnpye 1158 METHODS
213     (* 'specify' inherited from propagation is correct *)
214 aw0a 1 END trapezoid;
215    
216 johnpye 1158 (*-------------------------------
217     Midpoint method
218     *)
219 aw0a 1 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 johnpye 1158 (* 'specify' inherited from propagation is correct *)
244 aw0a 1 END midpoint;
245    
246 johnpye 1158 (*-------------------------------
247     Simpson's method
248     *)
249 aw0a 1 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 johnpye 1158 (* 'specify' inherited from propagation is correct *)
278 aw0a 1 END simpsons;
279    
280 johnpye 1158 (*-------------------------------
281     4-point Runge-Kutta-Gill method
282     *)
283 aw0a 1 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 johnpye 1158 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 aw0a 1
317     METHODS
318 johnpye 1158 (* 'specify' inherited from propagation is correct *)
319 aw0a 1 END runge_kutta_gill;
320    
321 johnpye 1158 (*-------------------------------
322     4-point Gauss method
323     *)
324 aw0a 1 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 johnpye 1158 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 aw0a 1
361     METHODS
362 johnpye 1158 (* 'specify' inherited from propagation is correct *)
363 aw0a 1 END gauss;
364    
365 johnpye 1158 (*-------------------------------
366     Lobatto method
367     *)
368 aw0a 1 MODEL lobatto(
369 johnpye 1158 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 aw0a 1
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 johnpye 1158 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 aw0a 1
411     METHODS
412 johnpye 1158 (* 'specify' inherited from propagation is correct *)
413 aw0a 1 END lobatto;
414    
415 johnpye 1158 (*-------------------------------
416     INTEGRATION model
417     *)
418 aw0a 1 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 johnpye 1158 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 aw0a 1 END integration;
522    
523 johnpye 1158 (*-------------------------------
524     Demo and test models
525     *)
526    
527 aw0a 1 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 johnpye 1158 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 aw0a 1 END test_dynamics;
548    
549 johnpye 1158 (*-------------------------------
550     'dynamic point'
551     *)
552 aw0a 1 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 johnpye 1158 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 aw0a 1 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 johnpye 1158 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 johnpye 1198 METHOD on_load;
617     RUN reset; RUN default_self; RUN values;
618     END on_load;
619 johnpye 1261 METHOD self_test;
620     (* no tests yet *)
621     END self_test;
622 aw0a 1 END bvp_test;
623 johnpye 1158 (* :ex: set ts=4: *)

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