/[ascend]/trunk/models/test/transamp.a4c
ViewVC logotype

Contents of /trunk/models/test/transamp.a4c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1254 - (show annotations) (download) (as text)
Sat Jan 27 10:50:40 2007 UTC (17 years, 5 months ago) by johnpye
File MIME type: text/x-ascend
File size: 3147 byte(s)
Added transamp.a4c, another one of the tests from Univ. of Bari.
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 MODEL transamp;
18
19 U_b, U_F, alpha, beta, C[1..5], R[0..9] IS_A real_constant;
20 FOR k IN [1..5] CREATE
21 C[k] :== k * 1e-6;
22 END FOR;
23
24 R[0] :== 1000;
25 FOR k IN [1..9] CREATE
26 R[k] :== 9000;
27 END FOR;
28
29 U_b :== 6;
30 U_F :== 0.026;
31 alpha :== 0.99;
32 beta :== 1e-6;
33
34 y[1..8], dy_dt[1..8] IS_A solver_var;
35 t IS_A solver_var;
36
37 U_e IS_A solver_var;
38 U_e = 0.1 * sin(200{PI}*t);
39
40 g23, g56 IS_A solver_var;
41 (* note, comments about possible overflow with this fn *)
42 g23 = beta * (exp((y[2] - y[3])/U_F) - 1);
43 g56 = beta * (exp((y[5] - y[6])/U_F) - 1);
44
45 -C[1]*dy_dt[1] + C[1]*dy_dt[2] = (y[1] - U_e) / R[0];
46 C[1]*dy_dt[1] - C[1]*dy_dt[2] = y[2] / R[1] + (y[2] - U_b) / R[2] + (1 - alpha)*g23;
47 -C[2]*dy_dt[3] = y[3] / R[3] - g23;
48 -C[3]*dy_dt[4] + C[3]*dy_dt[5] = (y[4] - U_b) / R[4] + alpha*g23;
49 C[3]*dy_dt[4] - C[3]*dy_dt[5] = y[5] / R[5] + (y[5] - U_b) / R[6] + (1 - alpha)*g56;
50 -C[4]*dy_dt[6] = y[6] / R[7] - g56;
51 -C[5]*dy_dt[7] + C[5]*dy_dt[8] = (y[7] - U_b) / R[8] + alpha*g56;
52 C[5]*dy_dt[7] - C[5]*dy_dt[8] = y[8] / R[9];
53
54 METHODS
55 METHOD values;
56 y[1] := 0;
57 y[2] := U_b/(R[2]/R[1] + 1);
58 y[3] := U_b/(R[2]/R[1] + 1);
59 y[4] := U_b;
60 y[5] := U_b/(R[6]/R[5] + 1);
61 y[6] := U_b/(R[6]/R[5] + 1);
62 y[7] := U_b;
63 y[8] := 0;
64
65 dy_dt[3] := -y[2]/(C[2]*R[3]);
66 dy_dt[6] := -y[5]/(C[4]*R[7]);
67
68 (* these were already calculated, numerically *)
69 dy_dt[1] := 51.338775;
70 dy_dt[2] := 51.338775;
71 dy_dt[4] := -24.9757667;
72 dy_dt[5] := -24.9757667;
73 dy_dt[7] := -10.00564453;
74 dy_dt[8] := -10.00564453;
75
76 t := 0 {s};
77 END values;
78
79 METHOD specify;
80 FIX t; (* t is incident in the equations, so we need to explicitly fix it *)
81 END specify;
82
83 METHOD ode_init;
84 FOR i IN [1..8] DO
85 y[i].ode_id := i; y[i].ode_type := 1;
86 dy_dt[i].ode_id := i; dy_dt[i].ode_type := 2;
87 END FOR;
88 FOR i IN [1..8] DO
89 y[i].obs_id := i;
90 END FOR;
91 t.ode_type := -1;
92 END ode_init;
93
94 METHOD on_load;
95 RUN reset;
96 RUN values;
97 RUN ode_init;
98 END on_load;
99
100 METHOD self_test;
101 ASSERT abs( y[1] + 0.5562145012262709e-002 ) < 1e-7;
102 ASSERT abs( y[2] - 3.006522471903042 ) < 1e-7;
103 ASSERT abs( y[3] - 2.849958788608128 ) < 1e-7;
104 ASSERT abs( y[4] - 2.926422536206241 ) < 1e-7;
105 ASSERT abs( y[5] - 2.704617865010554 ) < 1e-7;
106 ASSERT abs( y[6] - 2.761837778393145 ) < 1e-6;
107 ASSERT abs( y[7] - 4.770927631616772 ) < 1e-7;
108 ASSERT abs( y[8] - 1.236995868091548 ) < 1e-7;
109 END self_test;
110 END transamp;

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