Hi all,
[This is quite long, tl;dr: I cannot solve a matrix in reasonable time
that can be solved with Maxima in 8s.]
I'm using Lcapy (which in turn uses sympy) for computing a transfer
function of a simple circuit (see below for an ascii-art drawing).
I was referred here by the author of Lcapy.
As many of you may be aware one can compute the parameters of a circuit
by setting up a matrix and solving for the voltages (or currents). See
at the end for an ascii-art of the circuit.
Lcapy comes up with the following matrix:
from sympy import Matrix, symbols
R1, R2, C1, C2, C3, C4, C5, C6, C7, L1, L2, L3 = symbols(
'R1, R2, C1, C2, C3, C4, C5, C6, C7, L1, L2, L3', positive=True,
real=True)
s = symbols('s')
b = Matrix([0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1])
A = Matrix([
[1/R1, -1/R1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1],
[-1/R1, C1*s + 1/R1, -C1*s, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, -C1*s, C1*s + C2*s + C3*s, -C2*s, -C3*s, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, -C2*s, C2*s, 0, 0, 0, 0, 0, 1, 0, 0, 0],
[0, 0, -C3*s, 0, C3*s + C4*s + C5*s, -C4*s, -C5*s, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, -C4*s, C4*s, 0, 0, 0, 0, 1, 0, 0],
[0, 0, 0, 0, -C5*s, 0, C5*s + C6*s + C7*s, -C7*s, -C6*s, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, -C7*s, C7*s, 0, 0, 0, 1, 0],
[0, 0, 0, 0, 0, 0, -C6*s, 0, C6*s + 1/R2, 0, 0, 0, 0],
[0, 0, 0, 1, 0, 0, 0, 0, 0, -L1*s, 0, 0, 0],
[0, 0, 0, 0, 0, 1, 0, 0, 0, 0, -L2*s, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, -L3*s, 0],
[1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]])
x = A.solve(b, method='LU')
Trying to solve this with the 'LU' solver produces a really huge
solution (but the solution is quite fast):
def exsize (expr):
n = 0
if not expr.args:
return 1
for arg in expr.args:
n += exsize (arg)
return n
>>> print ([exsize (k) for k in x])
[52298, 26148, 884884, 128156, 730566, 65218, 510990, 35741, 255579,
128154, 65216, 35739, 26146]
[Is there a better way to determine the size of an expression?]
Simplifying the first expression (x [0]) should yield 1 (probably using
scipy.cancel). But I'm too impatient to wait for it, I've killed it when
it took 1.2GB resident memory after an hour or so.
Other solvers are *very* slow. Trying to solve with default 'GJ' method
wasn't completed in several hours, so I left it running overnight and
had terminated in the morning, unfortunately I had not computed the
sizes as above.
The size of the result created by 'LU' prevents further computations to
happen in reasonable time.
What I'm trying to achieve here is to compute the transfer function "H"
of the circuit, this is V(9)/V(1) (see circuit for points). It is x [8]
in the solution above (if it could be computed with sympy). I can do
this in sympy as below (without computing the voltages using the matrix,
instead the solution involves viewing the circuit as a cascade of
voltage dividers) in less than a second. The result is of reasonable
size (less than 1k).
I can successfully solve the matrix above in Maxima in 8 seconds
yielding a reasonably-sized result that can be further processed (with
ratsimp in my case). With ratsimp (x [1]) I indeed get 1 (indexes are
1-based in Maxima). And x [9] yields the correct transfer function.
Now my question: Should I expect such a huge matrix solution?
Is there another solver and/or another way to produce a reasonably-sized
result in reasonable time (several seconds)? Or can I force repeated
cancellation during matrix solving to prevent it from getting so big?
The state of the art in all intro texts on circuit simulation seems to
be to solve a matrix for the currents or voltages, should we use a
different method here (like the voltage divider solution below which is
applicable only to that particular ladder topology)?
Fast code to compute transfer function H for circuit below using voltage
dividers for the ladder components:
from sympy import solve, symbols
R1, R2, C1, C2, C3, C4, C5, C6, C7, L1, L2, L3, U1 = symbols(
'R1, R2, C1, C2, C3, C4, C5, C6, C7, L1, L2, L3, U1', positive=True,
real=True)
s = symbols('s')
R_V = R1 + 1 / (s * C1)
R_U7 = 1/(1/(1/(s*C7)+s*L3)+1/(1/(s*C6)+R2))
R_U5 = 1/(1/(1/(s*C4)+s*L2)+1/(1/(s*C5)+R_U7))
R_U3 = 1/(1/(1/(s*C2)+s*L1)+1/(1/(s*C3)+R_U5))
U3 = (R_U3 / (R_V + R_U3)) * U1
U5 = (R_U5 / (1 / (s * C3) + R_U5)) * U3
U7 = (R_U7 / (1 / (s * C5) + R_U7)) * U5
U9 = (R2 / (1 / (s * C6) + R2)) * U7
H = U9 / U1
[This should be identical to x [8] in the matrix solution above which so
far I've not been able to compute with sympy in a useable form]
Circuit ASCII-Art:
1 +----+ 2 C1 3 C3 5 C5 7 C6 9
o--| R1 |--+--| |---+---| |---+---| |---+---| |---+-----o
+----+ | | | |
| C2 | C4 | C7 |
--- --- --- ---
--- --- --- | | R2
|4 |6 |8 | |
> > > ---
> L1 > L2 > L3 |
> > > |
| | | |
o-------------------+---------+---------+---------+-----o
Thanks
Ralf
--
Dr. Ralf Schlatterbeck Tel: +43/2243/26465-16
Open Source Consulting www: www.runtux.com
Reichergasse 131, A-3411 Weidling email: [email protected]
--
You received this message because you are subscribed to the Google Groups
"sympy" group.
To unsubscribe from this group and stop receiving emails from it, send an email
to [email protected].
To view this discussion on the web visit
https://groups.google.com/d/msgid/sympy/20230825090051.oqpst544d6icdp4p%40runtux.com.
