models/test/transamp.a4c

REQUIRE "ivpsystem.a4l";
REQUIRE "atoms.a4l";
(*
	This is the 'Transistor amplifier' problem from the Test Set for IVP Solvers
	
	The problem is a stiff Differential Algebraic Equation of index 1, 
	consisting of 8 equations. The problem originates from electrical circuit 
	analysis, and is a model for the transistor amplifier. A complete 
	description of the problem could be found in transamp.pdf

	http://pitagora.dm.uniba.it/~testset/problems/transamp.php

	We're not really concerned with the details of this problem. We're just	
	trying to determine that our solver gets the right answer. High-precision
	solution results are provided online and we aim to reproduce those here.

	Note that, according to the graphs given at the above link, the DASSL
	solver is not so great at this particular problem. Not sure why this is,
	but it was certainly observed that it wasn't possible to obtain results
	at the same precision as obtained for the other examples from the Test Set.
*)
MODEL transamp;

	U_b, U_F, alpha, beta, C[1..5], R[0..9] IS_A real_constant;
	FOR k IN [1..5] CREATE
		C[k] :== k * 1e-6;
	END FOR;
	
	R[0] :== 1000;
	FOR k IN [1..9] CREATE
		R[k] :== 9000;
	END FOR;

	U_b :== 6;
	U_F :== 0.026;
	alpha :== 0.99;
	beta :== 1e-6;
	
	y[1..8], dy_dt[1..8] IS_A solver_var;
	t IS_A solver_var;

	U_e IS_A solver_var;
	U_e = 0.1 * sin(200{PI}*t);

	g23, g56 IS_A solver_var;
	(* note, comments about possible overflow with this fn *)
	g23 = beta * (exp((y[2] - y[3])/U_F) - 1);
	g56 = beta * (exp((y[5] - y[6])/U_F) - 1);

	-C[1]*dy_dt[1] + C[1]*dy_dt[2] = (y[1] - U_e) / R[0];
	 C[1]*dy_dt[1] - C[1]*dy_dt[2] = y[2] / R[1] + (y[2] - U_b) / R[2] + (1 - alpha)*g23;
	-C[2]*dy_dt[3]                 = y[3] / R[3] - g23;
	-C[3]*dy_dt[4] + C[3]*dy_dt[5] = (y[4] - U_b) / R[4] + alpha*g23;
	 C[3]*dy_dt[4] - C[3]*dy_dt[5] = y[5] / R[5] + (y[5] - U_b) / R[6] + (1 - alpha)*g56;
	-C[4]*dy_dt[6]                 = y[6] / R[7] - g56;
	-C[5]*dy_dt[7] + C[5]*dy_dt[8] = (y[7] - U_b) / R[8] + alpha*g56;
	 C[5]*dy_dt[7] - C[5]*dy_dt[8] = y[8] / R[9];

METHODS
	METHOD values;
		y[1] := 0;
		y[2] := U_b/(R[2]/R[1] + 1);
		y[3] := U_b/(R[2]/R[1] + 1);
		y[4] := U_b;
		y[5] := U_b/(R[6]/R[5] + 1);
		y[6] := U_b/(R[6]/R[5] + 1);
		y[7] := U_b;
		y[8] := 0;	

		dy_dt[3] := -y[2]/(C[2]*R[3]);
		dy_dt[6] := -y[5]/(C[4]*R[7]);

		(* these were already calculated, numerically *)
		dy_dt[1] :=  51.338775;
		dy_dt[2] :=  51.338775;
		dy_dt[4] := -24.9757667;
		dy_dt[5] := -24.9757667;
		dy_dt[7] := -10.00564453;
		dy_dt[8] := -10.00564453;

		t := 0 {s};
	END values;

	METHOD specify;
		FIX t; (* t is incident in the equations, so we need to explicitly fix it *)
	END specify;

	METHOD ode_init;
		FOR i IN [1..8] DO
			y[i].ode_id := i; y[i].ode_type := 1;
			dy_dt[i].ode_id := i; dy_dt[i].ode_type := 2;
		END FOR;
		FOR i IN [1..8] DO
			y[i].obs_id := i;
		END FOR;
		t.ode_type := -1;
	END ode_init;

	METHOD on_load;
		RUN reset;
		RUN values;
		RUN ode_init;
	END on_load;

	METHOD self_test;
		ASSERT abs( y[1] +  0.5562145012262709e-002 ) < 1e-7;
		ASSERT abs( y[2] -  3.006522471903042 ) < 1e-7;
		ASSERT abs( y[3] -  2.849958788608128 ) < 1e-7;
		ASSERT abs( y[4] -  2.926422536206241 ) < 1e-7;
		ASSERT abs( y[5] -  2.704617865010554 ) < 1e-7;
		ASSERT abs( y[6] -  2.761837778393145 ) < 1e-6;
		ASSERT abs( y[7] -  4.770927631616772 ) < 1e-7;
		ASSERT abs( y[8] -  1.236995868091548 ) < 1e-7;	
	END self_test;
END transamp;
1	REQUIRE "ivpsystem.a4l";
2	REQUIRE "atoms.a4l";
3	(*
4	This is the 'Transistor amplifier' problem from the Test Set for IVP Solvers
5
6	The problem is a stiff Differential Algebraic Equation of index 1,
7	consisting of 8 equations. The problem originates from electrical circuit
8	analysis, and is a model for the transistor amplifier. A complete
9	description of the problem could be found in transamp.pdf
10
11	http://pitagora.dm.uniba.it/~testset/problems/transamp.php
12
13	We're not really concerned with the details of this problem. We're just
14	trying to determine that our solver gets the right answer. High-precision
15	solution results are provided online and we aim to reproduce those here.
16
17	Note that, according to the graphs given at the above link, the DASSL
18	solver is not so great at this particular problem. Not sure why this is,
19	but it was certainly observed that it wasn't possible to obtain results
20	at the same precision as obtained for the other examples from the Test Set.
21	*)
22	MODEL transamp;
23
24	U_b, U_F, alpha, beta, C[1..5], R[0..9] IS_A real_constant;
25	FOR k IN [1..5] CREATE
26	C[k] :== k * 1e-6;
27	END FOR;
28
29	R[0] :== 1000;
30	FOR k IN [1..9] CREATE
31	R[k] :== 9000;
32	END FOR;
33
34	U_b :== 6;
35	U_F :== 0.026;
36	alpha :== 0.99;
37	beta :== 1e-6;
38
39	y[1..8], dy_dt[1..8] IS_A solver_var;
40	t IS_A solver_var;
41
42	U_e IS_A solver_var;
43	U_e = 0.1 * sin(200{PI}*t);
44
45	g23, g56 IS_A solver_var;
46	(* note, comments about possible overflow with this fn *)
47	g23 = beta * (exp((y[2] - y[3])/U_F) - 1);
48	g56 = beta * (exp((y[5] - y[6])/U_F) - 1);
49
50	-C[1]dy_dt[1] + C[1]dy_dt[2] = (y[1] - U_e) / R[0];
51	C[1]dy_dt[1] - C[1]dy_dt[2] = y[2] / R[1] + (y[2] - U_b) / R[2] + (1 - alpha)*g23;
52	-C[2]*dy_dt[3] = y[3] / R[3] - g23;
53	-C[3]dy_dt[4] + C[3]dy_dt[5] = (y[4] - U_b) / R[4] + alpha*g23;
54	C[3]dy_dt[4] - C[3]dy_dt[5] = y[5] / R[5] + (y[5] - U_b) / R[6] + (1 - alpha)*g56;
55	-C[4]*dy_dt[6] = y[6] / R[7] - g56;
56	-C[5]dy_dt[7] + C[5]dy_dt[8] = (y[7] - U_b) / R[8] + alpha*g56;
57	C[5]dy_dt[7] - C[5]dy_dt[8] = y[8] / R[9];
58
59	METHODS
60	METHOD values;
61	y[1] := 0;
62	y[2] := U_b/(R[2]/R[1] + 1);
63	y[3] := U_b/(R[2]/R[1] + 1);
64	y[4] := U_b;
65	y[5] := U_b/(R[6]/R[5] + 1);
66	y[6] := U_b/(R[6]/R[5] + 1);
67	y[7] := U_b;
68	y[8] := 0;
69
70	dy_dt[3] := -y[2]/(C[2]*R[3]);
71	dy_dt[6] := -y[5]/(C[4]*R[7]);
72
73	(* these were already calculated, numerically *)
74	dy_dt[1] := 51.338775;
75	dy_dt[2] := 51.338775;
76	dy_dt[4] := -24.9757667;
77	dy_dt[5] := -24.9757667;
78	dy_dt[7] := -10.00564453;
79	dy_dt[8] := -10.00564453;
80
81	t := 0 {s};
82	END values;
83
84	METHOD specify;
85	FIX t; (* t is incident in the equations, so we need to explicitly fix it *)
86	END specify;
87
88	METHOD ode_init;
89	FOR i IN [1..8] DO
90	y[i].ode_id := i; y[i].ode_type := 1;
91	dy_dt[i].ode_id := i; dy_dt[i].ode_type := 2;
92	END FOR;
93	FOR i IN [1..8] DO
94	y[i].obs_id := i;
95	END FOR;
96	t.ode_type := -1;
97	END ode_init;
98
99	METHOD on_load;
100	RUN reset;
101	RUN values;
102	RUN ode_init;
103	END on_load;
104
105	METHOD self_test;
106	ASSERT abs( y[1] + 0.5562145012262709e-002 ) < 1e-7;
107	ASSERT abs( y[2] - 3.006522471903042 ) < 1e-7;
108	ASSERT abs( y[3] - 2.849958788608128 ) < 1e-7;
109	ASSERT abs( y[4] - 2.926422536206241 ) < 1e-7;
110	ASSERT abs( y[5] - 2.704617865010554 ) < 1e-7;
111	ASSERT abs( y[6] - 2.761837778393145 ) < 1e-6;
112	ASSERT abs( y[7] - 4.770927631616772 ) < 1e-7;
113	ASSERT abs( y[8] - 1.236995868091548 ) < 1e-7;
114	END self_test;
115	END transamp;