Reputation: 2012
Following some example I found online, I can do this:
from sympy import var
from sympy import solve
Ldy, Ldz = var('Ldy Ldz')
g, x, y, z = var('g x y z')
xZ, yZ, zZ = var('xZ yZ zZ')
xdd, ydd, zdd = var('xdd ydd zdd')
E1 = z * xdd + (xZ - x) * (g + zdd)
E2 = z * ydd + (yZ - y) * (g + zdd) - Ldy
E3 = -y * xdd + x * ydd - zZ * (g + zdd) + Ldz
out = solve([E1, E2, E3], [xdd, ydd, Ldy])
print(type(xdd))
print("xdd = ", (out[xdd]).factor())
Which yields xdd = (g + zdd)*(x - xZ)/z
.
Now, doing it for my own equations:
from sympy import symbols, solve
x, y, z, k12, k26, x0 = symbols("x, y, z, k12, k26, x0")
symbols = x, y, z, k12, k26, x0
eq1 = k12 * x**2 -y
eq2 = k26 * y**3 - z
eq3 = x * 2*y + 6*z - x0
out = solve([eq1, eq2, eq3], [x,y,z])
print("x = ", (out[x]).factor())
Gives instead TypeError: list indices must be integers or slices, not Symbol
.
What am I doing wrong?
Upvotes: 1
Views: 132
Reputation:
The issue is that solve
has multiple return types: sometimes it returns a list, sometimes a dict, sometimes a list of dicts. The output form depends on the particulars of the equations being solved: number of variables, number of solutions. This mean one should use either list=True
or dict=True
to force consistent output from solve
. Note that dict=True
means the output is a list of dicts, since multiple solutions may exist -- which is the case here. In your example:
out = solve([eq1, eq2, eq3], [x,y,z], dict=True)
for sol in out:
print("x = ", sol[x].factor())
prints
x = 18**(1/3)*((3*x0 - sqrt(6*k12*k26*x0 + 1)/(k12*k26) + 1/(k12*k26))/k26)**(2/3)*(sqrt(6*k12*k26*x0 + 1) + 1)/(18*k12*x0)
x = -18**(1/3)*((3*x0 + sqrt(6*k12*k26*x0 + 1)/(k12*k26) + 1/(k12*k26))/k26)**(2/3)*(sqrt(6*k12*k26*x0 + 1) - 1)/(18*k12*x0)
x = -2**(1/3)*((3*x0 - sqrt(6*k12*k26*x0 + 1)/(k12*k26) + 1/(k12*k26))/k26)**(2/3)*(3**(2/3) - 3*3**(1/6)*I)*(sqrt(6*k12*k26*x0 +1) + 1)/(36*k12*x0)
x = -2**(1/3)*((3*x0 - sqrt(6*k12*k26*x0 + 1)/(k12*k26) + 1/(k12*k26))/k26)**(2/3)*(3**(2/3) + 3*3**(1/6)*I)*(sqrt(6*k12*k26*x0 +1) + 1)/(36*k12*x0)
x = 2**(1/3)*((3*x0 + sqrt(6*k12*k26*x0 + 1)/(k12*k26) + 1/(k12*k26))/k26)**(2/3)*(3**(2/3) - 3*3**(1/6)*I)*(sqrt(6*k12*k26*x0 + 1) - 1)/(36*k12*x0)
x = 2**(1/3)*((3*x0 + sqrt(6*k12*k26*x0 + 1)/(k12*k26) + 1/(k12*k26))/k26)**(2/3)*(3**(2/3) + 3*3**(1/6)*I)*(sqrt(6*k12*k26*x0 + 1) - 1)/(36*k12*x0)
For this and other reasons, SymPy developers recommend using solveset and its relatives instead of solve
. Specifically, nonlinsolve
can be used here:
out = nonlinsolve([eq1, eq2, eq3], [x,y,z])
for sol in out:
print("x = ", sol[x].factor())
which prints
x = -18**(1/3)*((3*x0 + sqrt(6*k12*k26*x0 + 1)/(k12*k26) + 1/(k12*k26))/k26)**(2/3)*(sqrt(6*k12*k26*x0 + 1) - 1)/(18*k12*x0)
x = 18**(1/3)*((3*x0 - sqrt(6*k12*k26*x0 + 1)/(k12*k26) + 1/(k12*k26))/k26)**(2/3)*(sqrt(6*k12*k26*x0 + 1) + 1)/(18*k12*x0)
x = 2**(1/3)*((3*x0 + sqrt(6*k12*k26*x0 + 1)/(k12*k26) + 1/(k12*k26))/k26)**(2/3)*(3**(2/3) - 3*3**(1/6)*I)*(sqrt(6*k12*k26*x0 + 1) - 1)/(36*k12*x0)
x = 2**(1/3)*((3*x0 + sqrt(6*k12*k26*x0 + 1)/(k12*k26) + 1/(k12*k26))/k26)**(2/3)*(3**(2/3) + 3*3**(1/6)*I)*(sqrt(6*k12*k26*x0 + 1) - 1)/(36*k12*x0)
x = -2**(1/3)*((3*x0 - sqrt(6*k12*k26*x0 + 1)/(k12*k26) + 1/(k12*k26))/k26)**(2/3)*(3**(2/3) + 3*3**(1/6)*I)*(sqrt(6*k12*k26*x0 +1) + 1)/(36*k12*x0)
x = -2**(1/3)*((3*x0 - sqrt(6*k12*k26*x0 + 1)/(k12*k26) + 1/(k12*k26))/k26)**(2/3)*(3**(2/3) - 3*3**(1/6)*I)*(sqrt(6*k12*k26*x0 +1) + 1)/(36*k12*x0)
The return type of solveset and its relatives is always a SymPy set.
Upvotes: 1