trunk/models/bvp.a4l

(*	ASCEND modelling environment
	Copyright (C) 1997 Benjamin A Allan
	Copyright (C) 1998, 2007 Carnegie Mellon University

	This program 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, or (at your option)
	any later version.

	This program is distributed in the 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 this program; if not, write to the Free Software
	Foundation, Inc., 59 Temple Place - Suite 330,
	Boston, MA 02111-1307, USA.
*)
REQUIRE "atoms.a4l";
(*
	This file contains models that implement an integration scheme completely	
	within ASCEND as a system of equations. Because data points are all 
	stored in memory, this approach uses a lot more memory (although with 	
	relation sharing, sparse matrix storage, etc, it's not completely terrible)

	@TODO more documentation here.

	@NOTE there is a test model implemented at the bottom of the file. See
	'test_dynamics' and 'bvp_test'. Note also that the model 'plotbvp.a4c'
	includes this test model and provides graphical output from the integration.

	@NODE some old comments were removed here as they confused the issue by
	referring to the LSODE integration engine.

	@NOTE The LSODE integration engine is provided as part of ASCEND and
	offers a more efficient way to solve integration models. Note that you will
	need to set up your model quite differently if you want to use that approach.

	Original code 1991.
	Decoupled from IVP code, rewritten in ASCEND IV -- Ben Allan, Mar 1997

	by Peter Piela, Boyd T. Safrit, Joseph J. Zaher, 1991
	remodelled by Benjamin A. Allan
*)

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.
	 	*)
		FIX x;
		FIX y[1..n_eq];
		FREE dydx[1..n_eq];
	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;
		FIX h;
		FREE final.x;
		FREE eval[1..n_point].y[1..n_eq];
		FREE eval[1..n_point].x;
	END specify;
END propagation;


(*-------------------------------
	Euler method
*)
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 *)
END euler;

(*-------------------------------
	Trapezoidal method
*)
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 *)
END trapezoid;

(*-------------------------------
	Midpoint method
*)
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 *)
END midpoint;

(*-------------------------------
	Simpson's method
*)
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 *)
END simpsons;

(*-------------------------------
	4-point Runge-Kutta-Gill method
*)
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 *)
END runge_kutta_gill;

(*-------------------------------
	4-point Gauss method
*)
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 *)
END gauss;

(*-------------------------------
	Lobatto method
*)
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 *)
END lobatto;

(*-------------------------------
	INTEGRATION model
*)
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;
		FIX stepsize;
		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;
		FIX x;
		FREE xdot;
		FIX t;
	END specify;
	METHOD values;
		xdot.lower_bound := -1e6{1/s};
	END values;
END test_dynamics;

(*-------------------------------
	'dynamic point'
*)
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;
		FIX stepsize;
		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;
(* :ex: set ts=4: *)
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.5h
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.5h
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.5h
305	eval[0].dydx[i];
306	eval[2].y[i] = eval[0].y[i] + 0.5h
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.11270167h
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.5h
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.88729833h
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.11270167h
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.5h
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.88729833h
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[(inpoint)..((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	END bvp_test;
619	(* :ex: set ts=4: *)