models/johnpye/idadenx.a4c

REQUIRE "ivpsystem.a4l";
REQUIRE "atoms.a4l";

(*
	ASCEND model based on the 'idadenx.c' example problem that accompanies IDA.
	The root-finding part isn't yet implemented though.
    ----------------------------------------------------------------------------
	This simple example problem for IDA, due to Robertson, 
	is from chemical kinetics, and consists of the following three 
	equations:

	     dy1/dt = -.04*y1 + 1.e4*y2*y3
	     dy2/dt = .04*y1 - 1.e4*y2*y3 - 3.e7*y2**2
	        0   = y1 + y2 + y3 - 1

	on the interval from t = 0.0 to t = 4.e10, with initial
	conditions: y1 = 1, y2 = y3 = 0.

	This model is tested (using the same parameters as used by idadenx.c) by the
	script pygtk/test.py.
*)
MODEL idadenx;
	y1, y2, y3 IS_A factor;
	dy1_dt IS_A factor;
	dy2_dt IS_A factor;

	eq1: dy1_dt = -0.04 * y1 + 1e4 * y2*y3;
	eq2: dy2_dt = +0.04 * y1 - 1e4 * y2*y3 - 3e7*y2^2;
	eq3: 0 = y1 + y2 + y3 - 1;
	
	t IS_A time;
METHODS

METHOD values; (* initial values *)
	y1 := 1;
	y2 := 0;
	y3 := 0;
	t := 0 {s};
END values;

METHOD specify;
	FIX y1, y2;
END specify;

METHOD ode_init;
	FREE y1, y2;
	t.ode_type := -1; t := 0 {s};

	dy1_dt.ode_id := 1; dy1_dt.ode_type := 2;
	y1.ode_id := 1; y1.ode_type := 1;

	dy2_dt.ode_id := 2; dy2_dt.ode_type := 2;
	y2.ode_id := 2; y2.ode_type := 1;

	(* y3.ode_id := 3; *)

	t.obs_id := 1;
	y1.obs_id := 2;
	y2.obs_id := 3;
	y3.obs_id := 4;
	dy1_dt.obs_id := 5;
	dy2_dt.obs_id := 6;

	y1.ode_atol := 1e-8;
	y2.ode_atol := 1e-14;
	y3.ode_atol := 1e-6;
END ode_init;
	
METHOD on_load;
	RUN reset; 
	RUN values;
	RUN ode_init;
END on_load;

END idadenx;
1	johnpye	942	REQUIRE "ivpsystem.a4l";
2			REQUIRE "atoms.a4l";
3
4	johnpye	951	(*
5			ASCEND model based on the 'idadenx.c' example problem that accompanies IDA.
6			The root-finding part isn't yet implemented though.
7			----------------------------------------------------------------------------
8			This simple example problem for IDA, due to Robertson,
9			is from chemical kinetics, and consists of the following three
10			equations:
11
12			dy1/dt = -.04y1 + 1.e4y2*y3
13			dy2/dt = .04y1 - 1.e4y2y3 - 3.e7y2**2
14			0 = y1 + y2 + y3 - 1
15
16			on the interval from t = 0.0 to t = 4.e10, with initial
17			conditions: y1 = 1, y2 = y3 = 0.
18
19			This model is tested (using the same parameters as used by idadenx.c) by the
20			script pygtk/test.py.
21			*)
22	johnpye	942	MODEL idadenx;
23			y1, y2, y3 IS_A factor;
24			dy1_dt IS_A factor;
25			dy2_dt IS_A factor;
26
27			eq1: dy1_dt = -0.04 * y1 + 1e4 * y2*y3;
28			eq2: dy2_dt = +0.04 * y1 - 1e4 * y2y3 - 3e7y2^2;
29			eq3: 0 = y1 + y2 + y3 - 1;
30
31			t IS_A time;
32			METHODS
33
34			METHOD values; (* initial values *)
35			y1 := 1;
36			y2 := 0;
37			y3 := 0;
38			t := 0 {s};
39			END values;
40
41			METHOD specify;
42			FIX y1, y2;
43			END specify;
44
45			METHOD ode_init;
46			FREE y1, y2;
47			t.ode_type := -1; t := 0 {s};
48	johnpye	1194
49	johnpye	942	dy1_dt.ode_id := 1; dy1_dt.ode_type := 2;
50			y1.ode_id := 1; y1.ode_type := 1;
51
52			dy2_dt.ode_id := 2; dy2_dt.ode_type := 2;
53			y2.ode_id := 2; y2.ode_type := 1;
54
55	johnpye	1193	(* y3.ode_id := 3; *)
56	johnpye	944
57			t.obs_id := 1;
58			y1.obs_id := 2;
59			y2.obs_id := 3;
60			y3.obs_id := 4;
61			dy1_dt.obs_id := 5;
62			dy2_dt.obs_id := 6;
63
64			y1.ode_atol := 1e-8;
65			y2.ode_atol := 1e-14;
66			y3.ode_atol := 1e-6;
67	johnpye	942	END ode_init;
68
69			METHOD on_load;
70	johnpye	1193	RUN reset;
71			RUN values;
72	johnpye	942	RUN ode_init;
73			END on_load;
74
75			END idadenx;