Commit | Line | Data |
---|---|---|
e890808f KB |
1 | (* |
2 | * Copyright (c) 1980 The Regents of the University of California. | |
3 | * All rights reserved. | |
4 | * | |
af359dea C |
5 | * Redistribution and use in source and binary forms, with or without |
6 | * modification, are permitted provided that the following conditions | |
7 | * are met: | |
8 | * 1. Redistributions of source code must retain the above copyright | |
9 | * notice, this list of conditions and the following disclaimer. | |
10 | * 2. Redistributions in binary form must reproduce the above copyright | |
11 | * notice, this list of conditions and the following disclaimer in the | |
12 | * documentation and/or other materials provided with the distribution. | |
13 | * 3. All advertising materials mentioning features or use of this software | |
14 | * must display the following acknowledgement: | |
15 | * This product includes software developed by the University of | |
16 | * California, Berkeley and its contributors. | |
17 | * 4. Neither the name of the University nor the names of its contributors | |
18 | * may be used to endorse or promote products derived from this software | |
19 | * without specific prior written permission. | |
e890808f | 20 | * |
af359dea C |
21 | * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND |
22 | * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE | |
23 | * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE | |
24 | * ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE | |
25 | * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL | |
26 | * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS | |
27 | * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) | |
28 | * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT | |
29 | * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY | |
30 | * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF | |
31 | * SUCH DAMAGE. | |
32 | * | |
33 | * @(#)parall.p 5.1 (Berkeley) 4/16/91 | |
e890808f KB |
34 | *) |
35 | ||
36 | program Parall(input,output); | |
37 | ||
38 | {Declare constants for unit conversions, convergence tests, etc.} | |
39 | ||
40 | const SQRT2 = 1.4142136; {Square root of 2} | |
41 | TWOPI = 6.2831853; {Two times pi} | |
42 | MINALPHA = 0.001; {Minimum y-step size} | |
43 | ROUGHLYZERO = 0.001; {Approximation to zero for convergence} | |
44 | YTHRESHOLD = 40.0; {Heuristic constant for thresholding} | |
45 | N = 8; {Vector and matrix dimension} | |
46 | ||
47 | ||
48 | {Declare types for circuit elements, vectors, and matrices.} | |
49 | ||
50 | type VSOURCE = record | |
51 | ampl: real; {Volts (peak)} | |
52 | freq: real; {Radians/second} | |
53 | xindex: integer; {Index for x value} | |
54 | yindex: integer; {Index for y value} | |
55 | end; | |
56 | ||
57 | RLPAIR = record | |
58 | r: real; {Ohms} | |
59 | l: real; {Henries} | |
60 | islope: real; {Amps/second} | |
61 | invariant: real; {Trapezoidal invariant} | |
62 | lasttime: real; {Most recent time} | |
63 | xindex: array [1..2] of integer; {x value indices} | |
64 | yindex: array [1..2] of integer; {y value indices} | |
65 | end; | |
66 | ||
67 | CAPACITOR = record | |
68 | c: real; {Farads} | |
69 | vslope: real; {Volts/second} | |
70 | invariant: real; {Trapezoidal invariant} | |
71 | lasttime: real; {Most recent time} | |
72 | xindex: array [1..2] of integer; {x value indices} | |
73 | yindex: array [1..2] of integer; {y value indices} | |
74 | end; | |
75 | ||
76 | VECTOR = array [1..N] of real; | |
77 | ||
78 | MATRIX = array [1..N,1..N] of real; | |
79 | ||
80 | ||
81 | {Declare global variables for central simulation information.} | |
82 | ||
83 | var ground: VSOURCE; {Ground -- a fake source} | |
84 | itcount: integer; {Main routine iteration count} | |
85 | update: integer; {Update loop counter for main} | |
86 | optimcount: integer; {Number of optimization steps} | |
87 | newtoncount: integer; {Number of Newton steps} | |
88 | v1: VSOURCE; {Voltage source V1} | |
89 | rl1: RLPAIR; {R1/L1 resistor-inductor pair} | |
90 | rl2: RLPAIR; {R2/L2 resistor-inductor pair} | |
91 | c1: CAPACITOR; {Capacitor C1} | |
92 | a: MATRIX; {Global matrix A} | |
93 | b: MATRIX; {Global matrix B} | |
94 | jac: MATRIX; {Global Jacobian matrix} | |
95 | x: VECTOR; {Global vector of dependents} | |
96 | xnext: VECTOR; {Next x-vector for simulation} | |
97 | y: VECTOR; {Global vector of independents} | |
98 | ynext: VECTOR; {Next y-vector for simulation} | |
99 | gradient: VECTOR; {Global gradient vector} | |
100 | h: real; {Time step value} | |
101 | time: real; {Current time} | |
102 | lastychange: real; {YStep's most recent y-change} | |
103 | timestep: integer; {Current timestep number} | |
104 | maxsteps: integer; {Number of time steps to run} | |
105 | oldxnorm: real; {Old one-norm of x-vector} | |
106 | newxnorm: real; {New one-norm of x-vector} | |
107 | closenough: boolean; {Flag to indicate convergence} | |
108 | ||
109 | ||
110 | ||
111 | ||
112 | {The following procedure initializes everything for the program based | |
113 | on the little test circuit suggested by Sarosh. The user is asked | |
114 | to specify the simulation and circuit parameters, and then the matrix | |
115 | and vector values are set up.} | |
116 | ||
117 | procedure InitializeEverything; | |
118 | var i,j: integer; | |
119 | rtemp: real; | |
120 | begin | |
121 | ||
122 | {Ready the input and output files (almost nil for Berkeley).} | |
123 | writeln(output); | |
124 | writeln(output,'*** Simulation Output Record ***'); | |
125 | writeln(output); | |
126 | writeln(output); | |
127 | ||
128 | {Initialize convergence test/indication variables.} | |
129 | oldxnorm := 0.0; | |
130 | newxnorm := 0.0; | |
131 | lastychange := 0.0; | |
132 | ||
133 | {Get desired time step size for simulation.} | |
134 | readln(input,h); | |
135 | writeln(output,'h (Seconds): ',h:12:7); | |
136 | ||
137 | {Get desired number of time steps for simulation.} | |
138 | readln(input,maxsteps); | |
139 | writeln(output,'maxsteps: ',maxsteps:4); | |
140 | ||
141 | {Get parameters for source V1 and initialize the source.} | |
142 | with v1 do | |
143 | begin | |
144 | readln(input,rtemp); | |
145 | writeln(output,'V1 (Volts RMS): ',rtemp:9:3); | |
146 | ampl := rtemp * SQRT2; | |
147 | readln(input,rtemp); | |
148 | writeln(output,'f (Hertz): ',rtemp:9:3); | |
149 | freq := rtemp * TWOPI; | |
150 | xindex := 1; | |
151 | yindex := 1; | |
152 | end; | |
153 | ||
154 | {Get parameters for R1/L1 pair and initialize the pair.} | |
155 | with rl1 do | |
156 | begin | |
157 | readln(input,r); | |
158 | writeln(output,'R1 (Ohms): ',r:9:3); | |
159 | readln(input,l); | |
160 | writeln(output,'L1 (Henries): ',l:12:7); | |
161 | islope := 0.0; | |
162 | invariant := 0.0; | |
163 | lasttime := -1.0; {Negative to force first update} | |
164 | xindex[1] := 2; | |
165 | yindex[1] := 2; | |
166 | xindex[2] := 3; | |
167 | yindex[2] := 3; | |
168 | end; | |
169 | ||
170 | {Get parameters for R2/L2 pair and initialize the pair.} | |
171 | with rl2 do | |
172 | begin | |
173 | readln(input,r); | |
174 | writeln(output,'R2 (Ohms): ',r:9:3); | |
175 | readln(input,l); | |
176 | writeln(output,'L2 (Henries): ',l:12:7); | |
177 | islope := 0.0; | |
178 | invariant := 0.0; | |
179 | lasttime := -1.0; {Negative to force first update} | |
180 | xindex[1] := 4; | |
181 | yindex[1] := 4; | |
182 | xindex[2] := 5; | |
183 | yindex[2] := 5; | |
184 | end; | |
185 | ||
186 | {Get parameters for capacitor C1 and initialize the device.} | |
187 | with c1 do | |
188 | begin | |
189 | readln(input,c); | |
190 | writeln(output,'C1 (Farads): ',c:12:7); | |
191 | vslope := 0.0; | |
192 | invariant := 0.0; | |
193 | lasttime := -1.0; {Negative to force first update} | |
194 | xindex[1] := 6; | |
195 | yindex[1] := 6; | |
196 | xindex[2] := 7; | |
197 | yindex[2] := 7; | |
198 | end; | |
199 | ||
200 | {Initialize the ground node.} | |
201 | with ground do | |
202 | begin | |
203 | ampl := 0.0; | |
204 | freq := 0.0; | |
205 | xindex := 8; | |
206 | yindex := 8; | |
207 | end; | |
208 | ||
209 | {Zero out all the vectors and matrices.} | |
210 | for i := 1 to N do | |
211 | begin | |
212 | x[i] := 0.0; | |
213 | y[i] := 0.0; | |
214 | for j := 1 to N do | |
215 | begin | |
216 | a[i,j] := 0.0; | |
217 | b[i,j] := 0.0; | |
218 | jac[i,j] := 0.0; | |
219 | end; | |
220 | end; | |
221 | ||
222 | {Initialize the A matrix.} | |
223 | a[1,2] := -1.0; | |
224 | a[2,3] := 1.0; | |
225 | a[2,4] := -1.0; | |
226 | a[3,5] := 1.0; | |
227 | a[4,7] := 1.0; | |
228 | a[5,1] := 1.0; | |
229 | a[7,6] := 1.0; | |
230 | a[8,8] := 1.0; | |
231 | ||
232 | {Initialize the B matrix.} | |
233 | b[1,1] := 1.0; | |
234 | b[3,7] := -1.0; | |
235 | b[4,8] := -1.0; | |
236 | b[5,2] := 1.0; | |
237 | b[6,3] := 1.0; | |
238 | b[6,4] := 1.0; | |
239 | b[7,5] := 1.0; | |
240 | b[8,6] := 1.0; | |
241 | ||
242 | {Initialize the Jacobian matrix.} | |
243 | rtemp := h / (2.0 * rl1.l + rl1.r * h); | |
244 | jac[2,2] := rtemp; | |
245 | jac[3,3] := rtemp; | |
246 | jac[2,3] := -rtemp; | |
247 | jac[3,2] := -rtemp; | |
248 | rtemp := h / (2.0 * rl2.l + rl2.r * h); | |
249 | jac[4,4] := rtemp; | |
250 | jac[5,5] := rtemp; | |
251 | jac[4,5] := -rtemp; | |
252 | jac[5,4] := -rtemp; | |
253 | jac[6,6] := -1.0; | |
254 | jac[7,7] := 1.0; | |
255 | jac[7,6] := h / (2.0 * c1.c); | |
256 | end; | |
257 | ||
258 | ||
259 | ||
260 | ||
261 | {The following procedure solves the equation Ax=b for an N x N system | |
262 | of linear equations, where A is the coefficient matrix, b is the | |
263 | right-hand-side vector, and x is the vector of unknowns. Gaussian | |
264 | elimination with maximal column pivots is used. } | |
265 | ||
266 | procedure Solve(a: MATRIX; b: VECTOR; var x: VECTOR); | |
267 | var y,z: real; | |
268 | i,j,k,k1: integer; | |
269 | begin | |
270 | for k := 1 to N-1 do | |
271 | begin | |
272 | y := abs(a[k,k]); | |
273 | j := k; | |
274 | k1 := k + 1; | |
275 | for i := k1 to N do | |
276 | if abs(a[i,k]) > y then | |
277 | begin | |
278 | j := i; | |
279 | y := abs(a[i,k]); | |
280 | end; | |
281 | for i := k to N do | |
282 | begin | |
283 | y := a[k,i]; | |
284 | a[k,i] := a[j,i]; | |
285 | a[j,i] := y; | |
286 | end; | |
287 | y := b[k]; | |
288 | b[k] := b[j]; | |
289 | b[j] := y; | |
290 | z := a[k,k]; | |
291 | for i := k1 to N do | |
292 | begin | |
293 | y := a[i,k] / z; | |
294 | a[i,k] := y; | |
295 | for j := k1 to N do a[i,j] := a[i,j] - y * a[k,j]; | |
296 | b[i] := b[i] - y * b[k]; | |
297 | end; | |
298 | end; | |
299 | x[N] := b[N] / a[N,N]; | |
300 | for i := N-1 downto 1 do | |
301 | begin | |
302 | y := b[i]; | |
303 | for j := i+1 to N do y := y - a[i,j] * x[j]; | |
304 | x[i] := y / a[i,i]; | |
305 | end; | |
306 | end; | |
307 | ||
308 | ||
309 | {The following procedure computes and returns a vector called "deltay", | |
310 | which is the change in the y-vector prescribed by the Newton-Rhapson | |
311 | algorithm.} | |
312 | ||
313 | procedure NewtonStep(var deltay: VECTOR); | |
314 | var phi: VECTOR; | |
315 | m: MATRIX; | |
316 | i,j,k: integer; | |
317 | begin | |
318 | for i := 1 to N do | |
319 | begin | |
320 | phi[i] := 0.0; | |
321 | for j := 1 to N do | |
322 | begin | |
323 | phi[i] := phi[i] + a[i,j] * y[j] + b[i,j] * x[j]; | |
324 | m[i,j] := -a[i,j]; | |
325 | for k := 1 to N do m[i,j] := m[i,j] - b[i,k] * jac[k,j]; | |
326 | end; | |
327 | ||
328 | end; | |
329 | Solve(m,phi,deltay); | |
330 | end; | |
331 | ||
332 | ||
333 | ||
334 | ||
335 | {The following function computes the value of theta, the objective | |
336 | function, given the x and y vectors.} | |
337 | ||
338 | function ThetaValue(x,y: VECTOR): real; | |
339 | var i,j: integer; | |
340 | phielem: real; | |
341 | theta: real; | |
342 | begin | |
343 | theta := 0.0; | |
344 | for i:= 1 to N do | |
345 | begin | |
346 | phielem := 0.0; | |
347 | for j := 1 to N do | |
348 | phielem := phielem + a[i,j] * y[j] + b[i,j] * x[j]; | |
349 | theta := theta + phielem * phielem; | |
350 | end; | |
351 | ThetaValue := theta; | |
352 | end; | |
353 | ||
354 | ||
355 | {The following function computes the theta value associated with a | |
356 | proposed step of size alpha in the direction of the gradient.} | |
357 | ||
358 | function Theta(alpha: real): real; | |
359 | var ythere: VECTOR; | |
360 | i: integer; | |
361 | begin | |
362 | for i := 1 to N do | |
363 | ythere[i] := y[i] - alpha * gradient[i]; | |
364 | Theta := ThetaValue(x,ythere); | |
365 | end; | |
366 | ||
367 | ||
368 | {The following procedure computes the gradient of the objective | |
369 | function (theta) with respect to the vector y.} | |
370 | ||
371 | procedure ComputeGradient; | |
372 | var i,j,k: integer; | |
373 | m: MATRIX; | |
374 | v: VECTOR; | |
375 | begin | |
376 | {Compute v = Ay + Bx and M = A' + J'B'.} | |
377 | for i := 1 to N do | |
378 | begin | |
379 | v[i] := 0.0; | |
380 | for j := 1 to N do | |
381 | begin | |
382 | v[i] := v[i] + a[i,j] * y[j] + b[i,j] * x[j]; | |
383 | m[i,j] := a[j,i]; | |
384 | for k := 1 to N do | |
385 | m[i,j] := m[i,j] + jac[k,i] * b[j,k]; | |
386 | end; | |
387 | end; | |
388 | {Compute gradient = 2Mv.} | |
389 | for i := 1 to N do | |
390 | begin | |
391 | gradient[i] := 0.0; | |
392 | for j := 1 to N do | |
393 | gradient[i] := gradient[i] + m[i,j] * v[j]; | |
394 | gradient[i] := 2.0 * gradient[i]; | |
395 | end; | |
396 | end; | |
397 | ||
398 | ||
399 | {The following procedure computes the bounds on alpha, the step size | |
400 | to take in the gradient direction. The bounding is done according | |
401 | to an algorithm suggested in S.W.Director's text on circuits.} | |
402 | ||
403 | procedure GetAlphaBounds(var lower,upper: real); | |
404 | var alpha: real; | |
405 | oldtheta,newtheta: real; | |
406 | begin | |
407 | if ThetaValue(x,y) = 0.0 | |
408 | then | |
409 | begin | |
410 | lower := 0; | |
411 | ||
412 | upper := 0; | |
413 | end | |
414 | else | |
415 | begin | |
416 | lower := MINALPHA; | |
417 | oldtheta := Theta(lower); | |
418 | upper := MINALPHA * 2.0; | |
419 | newtheta := Theta(upper); | |
420 | if newtheta <= oldtheta then | |
421 | begin | |
422 | alpha := upper; | |
423 | repeat | |
424 | begin | |
425 | oldtheta := newtheta; | |
426 | alpha := alpha * 2.0; | |
427 | newtheta := Theta(alpha); | |
428 | end | |
429 | until (newtheta > oldtheta); | |
430 | lower := alpha / 4.0; | |
431 | upper := alpha; | |
432 | end; | |
433 | end; | |
434 | end; | |
435 | ||
436 | ||
437 | {The following function computes the best value of alpha for the step | |
438 | in the gradient direction. This best value is the value that minimizes | |
439 | theta along the gradient-directed path.} | |
440 | ||
441 | function BestAlpha(lower,upper: real): real; | |
442 | var alphaL,alphaU,alpha1,alpha2: real; | |
443 | thetaL,thetaU,theta1,theta2: real; | |
444 | ||
445 | begin | |
446 | alphaL := lower; | |
447 | thetaL := Theta(alphaL); | |
448 | alphaU := upper; | |
449 | thetaU := Theta(alphaU); | |
450 | alpha1 := 0.381966 * alphaU + 0.618034 * alphaL; | |
451 | theta1 := Theta(alpha1); | |
452 | alpha2 := 0.618034 * alphaU + 0.381966 * alphaL; | |
453 | theta2 := Theta(alpha2); | |
454 | repeat | |
455 | if theta1 < theta2 | |
456 | then | |
457 | begin | |
458 | alphaU := alpha2; | |
459 | thetaU := theta2; | |
460 | alpha2 := alpha1; | |
461 | theta2 := theta1; | |
462 | alpha1 := 0.381966 * alphaU + 0.618034 * alphaL; | |
463 | theta1 := Theta(alpha1); | |
464 | end | |
465 | else | |
466 | begin | |
467 | alphaL := alpha1; | |
468 | thetaL := theta1; | |
469 | alpha1 := alpha2; | |
470 | theta1 := theta2; | |
471 | alpha2 := 0.618034 * alphaU + 0.381966 * alphaL; | |
472 | theta2 := Theta(alpha2); | |
473 | end | |
474 | until abs(thetaU - thetaL) <= ROUGHLYZERO; | |
475 | BestAlpha := (alphaL + alphaU) / 2.0; | |
476 | end; | |
477 | ||
478 | ||
479 | {The following procedure computes and returns a vector called "deltay", | |
480 | which is the change in the y-vector prescribed by the steepest-descent | |
481 | approach.} | |
482 | ||
483 | procedure OptimizationStep(var deltay: VECTOR); | |
484 | var lower,upper: real; | |
485 | alpha: real; | |
486 | i: integer; | |
487 | begin | |
488 | ComputeGradient; | |
489 | GetAlphaBounds(lower,upper); | |
490 | if lower <> upper then | |
491 | begin | |
492 | alpha := BestAlpha(lower,upper); | |
493 | for i:= 1 to N do deltay[i] := - alpha * gradient[i]; | |
494 | end; | |
495 | end; | |
496 | ||
497 | ||
498 | ||
499 | ||
500 | {The following function computes the one-norm of a vector argument. | |
501 | The length of the argument vector is assumed to be N.} | |
502 | ||
503 | function OneNorm(vec: VECTOR): real; | |
504 | var sum: real; | |
505 | i: integer; | |
506 | begin | |
507 | sum := 0; | |
508 | for i := 1 to N do sum := sum + abs(vec[i]); | |
509 | OneNorm := sum; | |
510 | end; | |
511 | ||
512 | ||
513 | {The following procedure takes a y-step, using the optimization | |
514 | approach when far from the solution and the Newton-Rhapson approach | |
515 | when fairly close to the solution.} | |
516 | ||
517 | procedure YStep; | |
518 | var deltay: VECTOR; | |
519 | ychange: real; | |
520 | scale: real; | |
521 | i: integer; | |
522 | begin | |
523 | NewtonStep(deltay); | |
524 | ychange := OneNorm(deltay); | |
525 | if ychange > YTHRESHOLD | |
526 | then | |
527 | { | |
528 | begin | |
529 | OptimizationStep(deltay); | |
530 | ychange := OneNorm(deltay); | |
531 | if ychange > YTHRESHOLD then | |
532 | } | |
533 | begin | |
534 | scale := YTHRESHOLD/ychange; | |
535 | for i := 1 to N do deltay[i] := scale * deltay[i]; | |
536 | optimcount := optimcount + 1; | |
537 | end {;} | |
538 | { | |
539 | optimcount := optimcount + 1; | |
540 | end | |
541 | } | |
542 | else | |
543 | begin | |
544 | newtoncount := newtoncount + 1; | |
545 | end; | |
546 | for i := 1 to N do ynext[i] := y[i] + deltay[i]; | |
547 | end; | |
548 | ||
549 | ||
550 | ||
551 | ||
552 | {The following procedure updates the output of a voltage source | |
553 | given the current time.} | |
554 | ||
555 | procedure VsourceStep(vn: VSOURCE); | |
556 | begin | |
557 | with vn do | |
558 | xnext[xindex] := ampl * sin(freq * time); | |
559 | end; | |
560 | ||
561 | ||
562 | {The following procedure updates the outputs of a resistor-inductor | |
563 | pair given the time step to take...that is, this routine takes a | |
564 | time step for resistor-inductor pairs. The new outputs are found | |
565 | by applying the trapezoidal rule.} | |
566 | ||
567 | procedure RLPairStep(var rln: RLPAIR); | |
568 | begin | |
569 | with rln do | |
570 | begin | |
571 | if (time > lasttime) then | |
572 | begin | |
573 | lasttime := time; | |
574 | invariant := xnext[xindex[1]] + (h / 2.0) * islope; | |
575 | end; | |
576 | islope := (y[yindex[1]] - y[yindex[2]] - r * xnext[xindex[1]]) / l; | |
577 | xnext[xindex[1]] := invariant + (h / 2.0) * islope; | |
578 | xnext[xindex[2]] := - xnext[xindex[1]]; | |
579 | end; | |
580 | end; | |
581 | ||
582 | ||
583 | {The following procedure updates the outputs of a capacitor given the | |
584 | time step...it takes the time step using a trapezoidal rule iteration.} | |
585 | ||
586 | procedure CapacitorStep(var cn: CAPACITOR); | |
587 | var v: real; | |
588 | begin | |
589 | with cn do | |
590 | begin | |
591 | if (time > lasttime) then | |
592 | begin | |
593 | lasttime := time; | |
594 | v := xnext[xindex[2]] - y[yindex[2]]; | |
595 | invariant := v + (h / 2.0) * vslope; | |
596 | end; | |
597 | vslope := y[yindex[1]] / c; | |
598 | v := invariant + (h / 2.0) * vslope; | |
599 | xnext[xindex[1]] := - y[yindex[1]]; | |
600 | xnext[xindex[2]] := y[yindex[2]] + v; | |
601 | end; | |
602 | end; | |
603 | ||
604 | ||
605 | {The following procedure controls the overall x-step for the | |
606 | specific circuit to be simulated.} | |
607 | ||
608 | procedure XStep; | |
609 | begin | |
610 | VsourceStep(v1); | |
611 | RLPairStep(rl1); | |
612 | RLPairStep(rl2); | |
613 | CapacitorStep(c1); | |
614 | VsourceStep(ground); | |
615 | end; | |
616 | ||
617 | ||
618 | ||
619 | ||
620 | {The following procedure prints out the values of all the voltages and | |
621 | currents for the circuit at each time step.} | |
622 | ||
623 | procedure PrintReport; | |
624 | begin | |
625 | writeln(output); | |
626 | writeln(output); | |
627 | writeln(output,'REPORT at Time = ',time:8:5,' seconds'); | |
628 | writeln(output,'Number of iterations used: ',itcount:2); | |
629 | writeln(output,'Number of truncations: ',optimcount:2); | |
630 | writeln(output,'Number of Newton y-steps: ',newtoncount:2); | |
631 | writeln(output,'The voltages and currents are:'); | |
632 | writeln(output,' V1 = ',x[1]:12:7,' I1 = ',y[1]:12:7); | |
633 | writeln(output,' V2 = ',y[2]:12:7,' I2 = ',x[2]:12:7); | |
634 | writeln(output,' V3 = ',y[3]:12:7,' I3 = ',x[3]:12:7); | |
635 | writeln(output,' V4 = ',y[4]:12:7,' I4 = ',x[4]:12:7); | |
636 | writeln(output,' V5 = ',y[5]:12:7,' I5 = ',x[5]:12:7); | |
637 | writeln(output,' V6 = ',x[7]:12:7,' I6 = ',y[6]:12:7); | |
638 | writeln(output,' V7 = ',y[7]:12:7,' I7 = ',x[6]:12:7); | |
639 | writeln(output,' V8 = ',x[8]:12:7,' I8 = ',y[8]:12:7); | |
640 | end; | |
641 | ||
642 | ||
643 | ||
644 | ||
645 | {Finally, the main routine controls the whole simulation process.} | |
646 | ||
647 | begin | |
648 | InitializeEverything; | |
649 | PrintReport; | |
650 | for timestep := 1 to maxsteps do | |
651 | begin | |
652 | itcount := 0; | |
653 | optimcount := 0; | |
654 | newtoncount := 0; | |
655 | closenough := FALSE; | |
656 | time := h * timestep; | |
657 | repeat | |
658 | begin | |
659 | itcount := itcount + 1; | |
660 | YStep; | |
661 | XStep; | |
662 | for update := 1 to N do | |
663 | begin | |
664 | x[update] := xnext[update]; | |
665 | y[update] := ynext[update]; | |
666 | end; | |
667 | oldxnorm := newxnorm; | |
668 | newxnorm := OneNorm(x); | |
669 | if abs(newxnorm - oldxnorm) <= ROUGHLYZERO | |
670 | then closenough := TRUE; | |
671 | end; | |
672 | if itcount < 4 then closenough := FALSE; | |
673 | until (itcount = 25) or closenough; | |
674 | PrintReport; | |
675 | end; | |
676 | end. | |
677 |