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

Contents of /trunk/models/bvp.a4l

Parent Directory Parent Directory | Revision Log Revision Log


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

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