0

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.

Ectobius Rex
  • 1
  • 1

1 Answers1

2

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)
smichr
  • 16,948
  • 2
  • 27
  • 34
  • Thanks a lot, it works just perfect! I can't upvote yet. But as soon I reach reputation I will. – Ectobius Rex Aug 31 '23 at 00:57
  • You can mark it as the correct answer to your question, even as a new contributor. – Malcolm Aug 31 '23 at 01:02
  • The reason this works is because your equations are numerically ill-conditioned with some very large and some very small coefficients. The `linsolve` solver will use partial pivoting if the coefficients are (inexact) floats and that can minimise the error for the given inexact inputs. That is enough in this case but might not be in some other cases. Alternatively consider using exact numbers at least for the arguments to `exp` which are what generate the large/small numbers. – Oscar Benjamin Aug 31 '23 at 07:32