trunk/models/bvp.a4l

REQUIRE "atoms.a4l";
(* => atoms.a4l, measures.a4l, system.a4l, basemodel.a4l *)
PROVIDE "bvp.a4l";

(*
 *  bvp.a4l
 *  by Benjamin A Allan
 *  remodeled from integration.a4l
 *  by Peter Piela, Boyd T. Safrit, Joseph J. Zaher
 *  Part of the ASCEND Library
 *  $Date: 1998/06/17 18:49:40 $
 *  $Revision: 1.6 $
 *  $Author: mthomas $
 *  $Source: /afs/cs.cmu.edu/project/ascend/Repository/models/bvp.a4l,v $
 *
 *  This file is part of the ASCEND Modeling Library.
 *
 *  Copyright (C) 1997 Benjamin A Allan
 *  Copyright (C) 1998 Carnegie Mellon University
 *
 *  The ASCEND Modeling Library is free software; you can redistribute
 *  it and/or modify it under the terms of the GNU General Public
 *  License as published by the Free Software Foundation; either
 *  version 2 of the License, or (at your option) any later version.
 *
 *  The ASCEND Modeling Library is distributed in hope that it
 *  will be useful, but WITHOUT ANY WARRANTY; without even the implied
 *  warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
 *  See the GNU General Public License for more details.
 *
 *  You should have received a copy of the GNU General Public License
 *  along with the program; if not, write to the Free Software
 *  Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139 USA.  Check
 *  the file named COPYING.
 *)

(*============================================================================*

    B V P . A 4 L
    -------------

    AUTHORS:     Ben Allan
                 Peter Piela
                 Boyd T. Safrit
                 Joseph J. Zaher

    DATES:       07/91 - Original code.
                 08/94 - Added the observations to LSODE.
                 09/94 - Added ivp as base type for all initial value problem
                         solvers. baa.
                 10/94 - Added other integration methods to supplement lsode
                         and provide a means for dynamic optimization.
                 04/96 - Modified for the use of constants and function abs
		 03/97 - Decoupled from IVP, rewritten in ASCEND IV.

    CONTENTS:    The bvp_point model definition that is required
                 in the integration model.
                 For boundary value
                 problems and for dynamic optimization, a general integration
                 model is supplied which enables users to specify any of a
                 number of alternative integration techniques.  The user
                 should be aware, however, that using the general integration
                 model will substantially increase the size of a problem.
                 When possible, an IVP library is preferable to use.

*)

MODEL bvp_base REFINES cmumodel;
NOTES 'purpose' SELF {
This MODEL anchors the boundary value problem library.
This library allows  you to solve differential/algebraic
equation systems using standard collocation methods.
See example bvp_test.
}
END NOTES;
END bvp_base;
(* *********************************************************** *)
(* ******  bvp_point: the collocation user MODEL interface *** *)
(* *********************************************************** *)
MODEL bvp_point REFINES bvp_base;
NOTES 'purpose' SELF {
Generic definition of the differential equation system,
  Dy[i] = dydx[i], i=1..n_eq, D = d/dx.
Create a refinement of this type that includes a part which
is your DAE MODEL and merge your variables with x,y,dydx as
needed. See MODEL bvp_test for an example.
Your refinement should define all the standard methods
default_self, etc.
}
'usage' SELF {
Note that the integration will behave much better if the
typical values of your state variables are all about the
same size. Integrators can lie if you give them poorly
scaled problems. Dimensionless variables are generally best
to use if you can easily formulate your MODEL in those
terms.
}
END NOTES;

	n_eq "number of implicit ODE's" IS_A integer_constant;

	dydx[1..n_eq] "derivative variables",
	y[1..n_eq] "state variables",
	x "independent variable"
		IS_A solver_var;

METHODS

METHOD default_self;
	(* do not assign anything as these are merged with
	 * user's variables which are all assigned.
	 *)
END default_self;

METHOD specify;
	(* The state variables are set
 	* TRUE, while their derivatives are FALSE.
 	*)
	x.fixed := TRUE;
	y[1..n_eq].fixed := TRUE;
	dydx[1..n_eq].fixed := FALSE;
END specify;
END bvp_point;


(* ******************************************************************** *)
(* *********************  bvp internal models ************************* *)
(* ******************************************************************** *)

MODEL propagation(
	n_point WILL_BE integer_constant;
	eval[0..n_point] WILL_BE bvp_point;
	n_eq WILL_BE integer_constant;
) WHERE (
	n_eq == eval[0].n_eq;
) REFINES bvp_base;
(*      eval[0..n_point] ARE_ALIKE;
 * Eval[i] must be deeply alike, so use of ARE_ALIKE here is heuristic at best.
 * Formal type compatibility does not ensure a well defined BVP.
 * No longer an issue, since eval is coming in from outside where it
 * was constructed all at the same time.
 *)

NOTES
'purpose' SELF {
This MODEL does nothing by itself, and should never be created.
Only the refinements of this MODEL should be created.
If you need a new collocation formula, create it by refining
propagation and then update the SELECT statement in
the MODEL integration.
}
'usage' n_order {
n_order did nothing.
It was mainly a carrier of information for the user.
Specifically, it gives the order of the error of the collocation.
in terms of h.
this information should just be kept in NOTES instead, like so.
}
END NOTES;


	h IS_A solver_var;
	(* propagate the dimensionality of x to h, though the
        h, eval[0].x ARE_ALIKE;
	 * ASCEND environment would do the same thing for us automagically.
	 * We actually couldn't care less about the type of h, so
	 * long as it is dimensionally consistent with the x.
	 * or we could do it with an assignment:
	h := eval[0].x;
	 *)
        initial ALIASES  eval[0];
        final ALIASES eval[n_point];
        final.x = initial.x + h;

METHODS
METHOD specify;
	RUN eval[0..n_point].specify;
	h.fixed := TRUE;
	final.x.fixed := FALSE;
	eval[1..n_point].y[1..n_eq].fixed := FALSE;
	eval[1..n_point].x.fixed := FALSE;
END specify;
END propagation;

MODEL euler(
	n_point WILL_BE integer_constant;
	eval[0..n_point] WILL_BE bvp_point;
	n_eq WILL_BE integer_constant;
) WHERE (
	n_eq == eval[0].n_eq;
	n_point == 1;
) REFINES propagation();

NOTES
'order of error' SELF {
The error in the euler method is O(h).
}
END NOTES;

        FOR i IN [1..n_eq] CREATE
            eval[1].y[i] = eval[0].y[i] + h*eval[0].dydx[i];
        END FOR;

METHODS
(* specify inherited from propagation is correct *)
(*
METHOD specify;
	RUN eval[0..n_point].specify;
	h.fixed := TRUE;
	final.x.fixed := FALSE;
	eval[1..n_point].y[1..n_eq].fixed := FALSE;
	eval[1..n_point].x.fixed := FALSE;
END specify;
*)

END euler;

MODEL trapezoid(
	n_point WILL_BE integer_constant;
	eval[0..n_point] WILL_BE bvp_point;
	n_eq WILL_BE integer_constant;
) WHERE (
	n_eq == eval[0].n_eq;
	n_point == 1;
) REFINES propagation();

NOTES
'order of error' SELF {
The error in the trapezoidal rule is O(h^2).
}
END NOTES;

        FOR i IN [1..n_eq] CREATE
            eval[1].y[i] = eval[0].y[i] + 0.5*h*
                (eval[0].dydx[i] + eval[1].dydx[i]);
        END FOR;

    METHODS
(* specify inherited from propagation is correct *)
(*
        METHOD specify;
            RUN eval[0..n_point].specify;
            h.fixed := TRUE;
            final.x.fixed := FALSE;
            eval[1..n_point].y[1..n_eq].fixed := FALSE;
            eval[1..n_point].x.fixed := FALSE;
        END specify;
*)

END trapezoid;

MODEL midpoint(
	n_point WILL_BE integer_constant;
	eval[0..n_point] WILL_BE bvp_point;
	n_eq WILL_BE integer_constant;
) WHERE (
	n_eq == eval[0].n_eq;
	n_point == 2;
) REFINES propagation();

NOTES
'order of error' SELF {
The error in the midpoint rule is O(h^2).
}
END NOTES;

        eval[1].x = initial.x + (1.0/2.0)*h;
        FOR i IN [1..n_eq] CREATE
            eval[1].y[i] = eval[0].y[i] + (1.0/2.0)*h*
                eval[0].dydx[i];
            eval[2].y[i] = eval[0].y[i] + h*
                eval[1].dydx[i];
        END FOR;

METHODS
(* specify inherited from propagation is correct *)
(*
METHOD specify;
	RUN eval[0..n_point].specify;
	h.fixed := TRUE;
	final.x.fixed := FALSE;
	eval[1..n_point].y[1..n_eq].fixed := FALSE;
	eval[1..n_point].x.fixed := FALSE;
END specify;
*)
END midpoint;

MODEL simpsons(
	n_point WILL_BE integer_constant;
	eval[0..n_point] WILL_BE bvp_point;
	n_eq WILL_BE integer_constant;
) WHERE (
	n_eq == eval[0].n_eq;
	n_point == 2;
) REFINES propagation();

NOTES
'order of error' SELF {
The error in simpson's rule is O(h^4).
}
END NOTES;

        eval[1].x = initial.x + 0.5*h;
        FOR i IN [1..n_eq] CREATE
            eval[1].y[i] = eval[0].y[i] + 0.5*h*
                ((5.0/12.0)*eval[0].dydx[i] +
                (8.0/12.0)*eval[1].dydx[i] -
                (1.0/12.0)*eval[2].dydx[i]);
            eval[2].y[i] = eval[0].y[i] + h*
                ((1.0/6.0)*eval[0].dydx[i] +
                (4.0/6.0)*eval[1].dydx[i] +
                (1.0/6.0)*eval[2].dydx[i]);
        END FOR;

METHODS
(* specify inherited from propagation is correct *)
(*
METHOD specify;
	RUN eval[0..n_point].specify;
	h.fixed := TRUE;
	final.x.fixed := FALSE;
	eval[1..n_point].y[1..n_eq].fixed := FALSE;
	eval[1..n_point].x.fixed := FALSE;
END specify;
*)

END simpsons;

MODEL runge_kutta_gill(
	n_point WILL_BE integer_constant;
	eval[0..n_point] WILL_BE bvp_point;
	n_eq WILL_BE integer_constant;
) WHERE (
	n_eq == eval[0].n_eq;
	n_point == 4;
) REFINES propagation();

NOTES
'order of error' SELF {
The error in the 4 point Runge-Kutta-Gill method is O(h^4).
}
END NOTES;

        eval[1].x = eval[2].x;
        eval[3].x = eval[4].x;
        eval[1].x = initial.x + 0.5*h;
        FOR i IN [1..n_eq] CREATE
            eval[1].y[i] = eval[0].y[i] + 0.5*h*
                eval[0].dydx[i];
            eval[2].y[i] = eval[0].y[i] + 0.5*h*
                (0.41421356*eval[0].dydx[i] +
                0.58578644*eval[1].dydx[i]);
            eval[3].y[i] = eval[0].y[i] + h*
                (-0.707106781*eval[1].dydx[i] +
                1.707106781*eval[2].dydx[i]);
            eval[4].y[i] = eval[0].y[i] + h*
                (0.166666667*eval[0].dydx[i] +
                0.097631073*eval[1].dydx[i] +
                0.569035590*eval[2].dydx[i] +
                0.166666667*eval[3].dydx[i]);
        END FOR;

METHODS
(* specify inherited from propagation is correct *)
(*
METHOD specify;
	RUN eval[0..n_point].specify;
	h.fixed := TRUE;
	final.x.fixed := FALSE;
	eval[1..n_point].y[1..n_eq].fixed := FALSE;
	eval[1..n_point].x.fixed := FALSE;
END specify;
*)

END runge_kutta_gill;

MODEL gauss(
	n_point WILL_BE integer_constant;
	eval[0..n_point] WILL_BE bvp_point;
	n_eq WILL_BE integer_constant;
) WHERE (
	n_eq == eval[0].n_eq;
	n_point == 4;
) REFINES propagation();

NOTES
'order of error' SELF {
The error in 4 point gauss method is O(h^6).
}
END NOTES;

        eval[1].x = initial.x + 0.11270167*h;
        eval[2].x = initial.x + 0.5*h;
        eval[3].x = initial.x + 0.88729833*h;
        FOR i IN [1..n_eq] CREATE
            eval[1].y[i] = initial.y[i] + 0.11270167*h*
                (1.23236*eval[1].dydx[i] -
                0.31922*eval[2].dydx[i] +
                0.08686*eval[3].dydx[i]);
            eval[2].y[i] = initial.y[i] + 0.5*h*
                (0.60053*eval[1].dydx[i] +
                0.44444*eval[2].dydx[i] -
                0.04497*eval[3].dydx[i]);
            eval[3].y[i] = initial.y[i] + 0.88729833*h*
                (0.30203*eval[1].dydx[i] +
                0.54144*eval[2].dydx[i] +
                0.15653*eval[3].dydx[i]);
            eval[4].y[i] = initial.y[i] + h*
                (0.27778*eval[1].dydx[i] +
                0.44444*eval[2].dydx[i] +
                0.27778*eval[3].dydx[i]);
        END FOR;

METHODS
(* specify inherited from propagation is correct *)
(*
METHOD specify;
	RUN eval[0..n_point].specify;
	h.fixed := TRUE;
	final.x.fixed := FALSE;
	eval[1..n_point].y[1..n_eq].fixed := FALSE;
	eval[1..n_point].x.fixed := FALSE;
specify;
*)

END gauss;

MODEL lobatto(
	n_point WILL_BE integer_constant;
	eval[0..n_point] WILL_BE bvp_point;
	n_eq WILL_BE integer_constant;
) WHERE (
	n_eq == eval[0].n_eq;
	n_point == 4;
) REFINES propagation();

NOTES
'order of error' SELF {
The error in the 4 lobatto point method is O(h^8).
}
END NOTES;

        eval[1].x = initial.x + 0.11270167*h;
        eval[2].x = initial.x + 0.5*h;
        eval[3].x = initial.x + 0.88729833*h;
        FOR i IN [1..n_eq] CREATE
            eval[1].y[i] = initial.y[i] + 0.11270167*h*
                (0.428010*eval[0].dydx[i] +
                0.602339*eval[1].dydx[i] -
                0.044301*eval[2].dydx[i] +
                0.029587*eval[3].dydx[i] -
                0.015635*eval[4].dydx[i]);
            eval[2].y[i] = initial.y[i] + 0.5*h*
                (-0.06250*eval[0].dydx[i] +
                0.68122*eval[1].dydx[i] +
                0.44444*eval[2].dydx[i] -
                0.12566*eval[3].dydx[i] +
                0.06250*eval[4].dydx[i]);
            eval[3].y[i] = initial.y[i] + 0.88729833*h*
                (0.001986*eval[0].dydx[i] +
                0.309303*eval[1].dydx[i] +
                0.506524*eval[2].dydx[i] +
                0.236552*eval[3].dydx[i] -
                0.054365*eval[4].dydx[i]);
            eval[4].y[i] = initial.y[i] + h*
                (0.277778*eval[1].dydx[i] +
                0.444444*eval[2].dydx[i] +
                0.277778*eval[3].dydx[i]);
        END FOR;

METHODS
(* specify inherited from propagation is correct *)
(*
METHOD specify;
	RUN eval[0..n_point].specify;
	h.fixed := TRUE;
	final.x.fixed := FALSE;
	eval[1..n_point].y[1..n_eq].fixed := FALSE;
	eval[1..n_point].x.fixed := FALSE;
END specify;
*)

END lobatto;

MODEL integration(
	nstep IS_A integer_constant;
	npoint IS_A integer_constant;
	method IS_A symbol_constant;
	grid[0..nstep*npoint] WILL_BE bvp_point;
	stepsize WILL_BE solver_var;
	n_eq IS_A integer_constant;
) WHERE (
	((method == 'Euler') AND (npoint == 1)) OR
	((method == 'trapezoidal') AND (npoint == 1)) OR
	((method == 'midpoint') AND (npoint == 2)) OR
	((method == 'Simpsons') AND (npoint == 2)) OR
	((method == 'Runge-Kutta-Gill') AND (npoint == 4)) OR
	((method == 'Gauss') AND (npoint == 4)) OR
	((method == 'Lobatto') AND (npoint == 4));
) REFINES bvp_base;

NOTES 'usage' SELF {
This writes the equations for a mesh with nstep collocation
intervals. In each interval there npoint nodes.
The total mesh has (nstep * npoint + 1) nodes.

In order to use this MODEL, create a refinement of bvp_point
that contains the MODEL you want to integrate, then create
an array of them (grid) and pass it in along with the other options.
See the MODEL bvp_point for details of how to do this.
See the MODELs bvp_test and bvp_dynamics for an example.
}
END NOTES;
	(* cook up the little arrays for propagation out of the
	 * big array from the user.
	 *)
	FOR i IN [0..nstep-1] CREATE
		step_nodes[i+1][pset[i]] ALIASES
			(grid[(i*npoint)..((i+1)*npoint)])
		WHERE pset[i] IS_A set OF integer_constant WITH_VALUE
			(0..npoint);
	END FOR;

	(* merge all the sets. we need a better way of
	 * getting at the set defined for an array.
	 * we could get by without the merge.
	 *)
	pset[0..nstep-1] ARE_THE_SAME;

	(* create the generic step stuctures and refine to rule requested *)
	FOR i IN [1..nstep] CREATE
		step[i] IS_A propagation(npoint,step_nodes[i],n_eq);
	END FOR;

	(* refine to rule *)
	SELECT (method)
	CASE 'Euler':
		FOR i IN [1..nstep] CREATE
    			step[i] IS_REFINED_TO
			euler(npoint,step_nodes[i],n_eq);
		END FOR;
	CASE 'trapezoidal':
		FOR i IN [1..nstep] CREATE
    			step[i] IS_REFINED_TO
			trapezoid(npoint,step_nodes[i],n_eq);
		END FOR;
	CASE 'midpoint':
		FOR i IN [1..nstep] CREATE
    			step[i] IS_REFINED_TO
			midpoint(npoint,step_nodes[i],n_eq);
		END FOR;
	CASE 'Simpsons':
		FOR i IN [1..nstep] CREATE
    			step[i] IS_REFINED_TO
			simpsons(npoint,step_nodes[i],n_eq);
		END FOR;
	CASE 'Runge-Kutta-Gill':
		FOR i IN [1..nstep] CREATE
    			step[i] IS_REFINED_TO
			runge_kutta_gill(npoint,step_nodes[i],n_eq);
		END FOR;
	CASE 'Gauss':
		FOR i IN [1..nstep] CREATE
    			step[i] IS_REFINED_TO
			gauss(npoint,step_nodes[i],n_eq);
		END FOR;
	CASE 'Lobatto':
		FOR i IN [1..nstep] CREATE
    			step[i] IS_REFINED_TO
			lobatto(npoint,step_nodes[i],n_eq);
		END FOR;
	END SELECT;


METHODS
METHOD default_self;
	RUN values;
END default_self;
METHOD values;
	step[1..nstep].h := stepsize;
END values;
METHOD specify;
	stepsize.fixed := TRUE;
	FOR i IN [1..nstep] DECREASING DO
		RUN step[i].specify;
	END FOR;
END specify;

END integration;

(************************************************************************)
(************** demo and test models ************************************)
(************************************************************************)
MODEL testbvp_base REFINES testcmumodel;
END testbvp_base;

MODEL test_dynamics REFINES testbvp_base;
	x IS_A factor;
	t IS_A time;
	xdot IS_A frequency;
	xdot = 2{1/s} * x * exp(-2{1/s}*t) * cos(t*1{rad/s});
METHODS
METHOD default_self;
	RUN values;
END default_self;
METHOD specify;
	x.fixed := TRUE;
	xdot.fixed := FALSE;
	t.fixed := TRUE;
END specify;
METHOD values;
	xdot.lower_bound := -1e6{1/s};
END values;
END test_dynamics;

MODEL dynamic_point REFINES bvp_point;

	n_eq :== 1; (* you can have as many as you like *)

	(* replace the following line with an IS_A of your MODEL *)
	your_model IS_A test_dynamics;

	(* connect independent variables. If your MODEL is autonomous,
	 * you don't need to, but then you must call
	 * bvp_point::clear and bvp_point::specify in your methods
	 * to get x set properly.
	 *)
	x, your_model.t ARE_THE_SAME;

	(* connect your variables to the integration state variables *)
	y[1], your_model.x ARE_THE_SAME;

	(* connect your variables to the integration derivatives *)
	dydx[1], your_model.xdot ARE_THE_SAME;

METHODS

METHOD specify;

	RUN your_model.specify;
	(* If you write your_model correctly, you don't need
	   to run bvp_point::specify because you've wired all
	   the bvp variables into some of yours.
	   (Not TRUE if your MODEL is autonomous, though).
		RUN bvp_point::specify;
	*)

END specify;
END dynamic_point;


MODEL bvp_test REFINES testbvp_base;
	stepsize IS_A time;
	npoint, nstep IS_A integer_constant;
	npoint :== 1;
	nstep :== 100;
	nodes[0..nstep*npoint] IS_A dynamic_point;
	initial ALIASES nodes[0];
	final ALIASES nodes[nstep*npoint];
	integral IS_A
	integration(nstep,npoint,'Euler',nodes,stepsize,initial.n_eq);

	(* we don't tie the test MODEL to the ASCEND plot MODEL *)

METHODS
METHOD default_self;
	stepsize := 0.05{s};
	initial.x := 0{s};
	initial.y[1] := 1;
	RUN nodes[0..nstep*npoint].your_model.default_self;
	RUN integral.default_self;
END default_self;
METHOD specify;
	stepsize.fixed := TRUE;
	RUN integral.specify;
END specify;
METHOD values;
	stepsize := 0.05{s};
	initial.x := 0{s};
	initial.y[1] := 1;
	RUN nodes[0..nstep*npoint].your_model.values;
	RUN integral.values;
END values;

END bvp_test;
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	x.fixed := TRUE;
120	y[1..n_eq].fixed := TRUE;
121	dydx[1..n_eq].fixed := FALSE;
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	h.fixed := TRUE;
179	final.x.fixed := FALSE;
180	eval[1..n_point].y[1..n_eq].fixed := FALSE;
181	eval[1..n_point].x.fixed := FALSE;
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	h.fixed := TRUE;
210	final.x.fixed := FALSE;
211	eval[1..n_point].y[1..n_eq].fixed := FALSE;
212	eval[1..n_point].x.fixed := FALSE;
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.5h
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	h.fixed := TRUE;
244	final.x.fixed := FALSE;
245	eval[1..n_point].y[1..n_eq].fixed := FALSE;
246	eval[1..n_point].x.fixed := FALSE;
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	h.fixed := TRUE;
281	final.x.fixed := FALSE;
282	eval[1..n_point].y[1..n_eq].fixed := FALSE;
283	eval[1..n_point].x.fixed := FALSE;
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.5h
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	h.fixed := TRUE;
321	final.x.fixed := FALSE;
322	eval[1..n_point].y[1..n_eq].fixed := FALSE;
323	eval[1..n_point].x.fixed := FALSE;
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.5h
349	eval[0].dydx[i];
350	eval[2].y[i] = eval[0].y[i] + 0.5h
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	h.fixed := TRUE;
369	final.x.fixed := FALSE;
370	eval[1..n_point].y[1..n_eq].fixed := FALSE;
371	eval[1..n_point].x.fixed := FALSE;
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.11270167h
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.5h
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.88729833h
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	h.fixed := TRUE;
420	final.x.fixed := FALSE;
421	eval[1..n_point].y[1..n_eq].fixed := FALSE;
422	eval[1..n_point].x.fixed := FALSE;
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.11270167h
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.5h
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.88729833h
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	h.fixed := TRUE;
477	final.x.fixed := FALSE;
478	eval[1..n_point].y[1..n_eq].fixed := FALSE;
479	eval[1..n_point].x.fixed := FALSE;
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[(inpoint)..((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	stepsize.fixed := TRUE;
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	x.fixed := TRUE;
608	xdot.fixed := FALSE;
609	t.fixed := TRUE;
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	stepsize.fixed := TRUE;
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;