I have to solve a system of equations that satisfice some conditions. When I do it with solve from Sympy, I only get 2 out 8 values correct.
Here's my code.
import numpy as np
from sympy import symbols, exp, Eq, solve
gamma0, gamma1, gamma2, gamma3, gamma4, gamma5, gamma6, gamma7 = symbols('gamma0 gamma1 gamma2 gamma3 gamma4 gamma5 gamma6 gamma7')
gamma = [gamma0, gamma1, gamma2, gamma3, gamma4, gamma5, gamma6, gamma7]
x = symbols('x')
a = 20
b = 5
u =[
[0.388518*gamma0*exp(-2.012490*x) - 0.675612*gamma1*exp(-0.369901*x) + 0.336576*gamma2*exp(0.369901*x) - 0.103503*gamma3*exp(2.012490*x) + 20.168718,
0.163057*gamma4*exp(-2.421891*x) - 0.951601*gamma5*exp(-0.951784*x) + 0.086427*gamma6*exp(0.754538*x) - 0.045759*gamma7*exp(2.528248*x) + 2.394974],
[-0.900008*gamma0*exp(-2.012490*x) - 0.520883*gamma1*exp(-0.369901*x) + 0.398678*gamma2*exp(0.369901*x) - 0.168303*gamma3*exp(2.012490*x) + 20.034468,
-0.976653*gamma4*exp(-2.421891*x) - 0.238418*gamma5*exp(-0.951784*x) + 0.122769*gamma6*exp(0.754538*x) - 0.080425*gamma7*exp(2.528248*x) + 2.329339],
[-0.168303*gamma0*exp(-2.012490*x) - 0.398678*gamma1*exp(-0.369901*x) + 0.520883*gamma2*exp(0.369901*x) - 0.900008*gamma3*exp(2.012490*x) + 19.932869,
-0.095726*gamma4*exp(-2.421891*x) - 0.120587*gamma5*exp(-0.951784*x) + 0.326236*gamma6*exp(0.754538*x) - 0.963388*gamma7*exp(2.528248*x) + 2.504797],
[-0.103503*gamma0*exp(-2.012490*x) - 0.336576*gamma1*exp(-0.369901*x) + 0.675612*gamma2*exp(0.369901*x) + 0.388518*gamma3*exp(2.012490*x) + 19.892514,
-0.101960*gamma4*exp(-2.421891*x) - 0.151889*gamma5*exp(-0.951784*x) + 0.933288*gamma6*exp(0.754538*x) + 0.251633*gamma7*exp(2.528248*x) + 3.335936]
]
from sympy import Eq, solve
# Define the symbolic equations
eqns = [
Eq(u[0][0].subs(x, 0), 0),
Eq(u[1][0].subs(x, 0), 0),
Eq(u[0][0].subs(x, b), u[0][1].subs(x, b)),
Eq(u[1][0].subs(x, b), u[1][1].subs(x, b)),
Eq(u[2][0].subs(x, b), u[2][1].subs(x, b)),
Eq(u[3][0].subs(x, b), u[3][1].subs(x, b)),
Eq(u[2][1].subs(x, a), 0),
Eq(u[3][1].subs(x, a), 0)
]
# Solve the equations for the gamma variables
gamma_symbols = [gamma0, gamma1, gamma2, gamma3, gamma4, gamma5, gamma6, gamma7]
solutions = solve(eqns, gamma_symbols)
# Print the solutions
for symbol, value in solutions.items():
print(f"{symbol}: {value}")
The values I need to obtain are the following:
gamma0: 3.3451862507144963721760546900288
gamma1: 30.086443529744653795656405805539
gamma2: -3.3917670835986188504875542362514
gamma3: 0.00013878499955117735347592615338105
gamma4: -775727.3319427781475689057795114
gamma5: -941.09415775171256674990468706789
gamma6: -0.0000010944407599058234322536821517657
gamma7: 0.00000000000000000000013958849298737787410748516872778
However when I run my code I get these values:
gamma0: -2953067.95180723
gamma1: -1383474.12048191
gamma2: -13623.4036144577
gamma3: 1.03375979503954
gamma4: -8212497765.78320
gamma5: -15039124.0279307
gamma6: -0.00000109854435279209
gamma7: 1.40162266686094E-22
Where only gamma6 and gamma7 are accurate. I wanna know if maybe Python can't handle the system or if there's another good way to solve it.
Definitely use linsolve for this:
>>> from sympy import *
# your code without solve
>>> v = [gamma0, gamma1, gamma2, gamma3, gamma4, gamma5, gamma6, gamma7]
>>> for i in zip(v, list(linsolve(eqns, v))[0]): # only one soln
... print(i)
...
(gamma0, 3.34518073050084)
(gamma1, 30.0864601804563)
(gamma2, -3.39178380795305)
(gamma3, 0.000138785046404537)
(gamma4, -775727.547533190)
(gamma5, -941.091400437524)
(gamma6, -1.09445480076177e-6)
(gamma7, 1.39590082130360e-22)