Nsolve will not solve

Question

I am trying to create a surface plot based on temperatures. I need to feed in a hot and cold temperature to a function that solves a system of equations for our "z-axis" value. The function works fine until I set it to some variable. The system doesn't give completely solved when I set it to a variable. Below is an example of the errors I get:

SympifyError                              Traceback (most recent call last)
<ipython-input-12-828bf02f4398> in <module>
     49 cin = linspace(0,200,100)
     50 X, Y = meshgrid(hin,cin)
---> 51 Z = solver(X,Y)
     52 
     53 ax = axes(projection='3d')

<ipython-input-12-828bf02f4398> in solver(TH, TC)
     34     Tinfhin = TH +273.15
     35     Tinfcin = TC + 273.15
---> 36     sols = sy.nsolve(  (Eq(Qh,mdoth * cph * (Tinfhin - Tinfhout) ),
     37           Eq(Qh,nsh * hh * Ash * ((Tinfhin + Tinfhout)/2 - Th)),
     38           Eq(Qh,n * (alpha * II * Th - 0.5 * (II**2) * ree + (Ke * (Th-Tc)))),

D:\Users\sampl\Anaconda3\lib\site-packages\sympy\core\relational.py in __new__(cls, lhs, rhs, **options)
    389 
    390         lhs = _sympify(lhs)
--> 391         rhs = _sympify(rhs)
    392 
    393         evaluate = options.pop('evaluate', global_evaluate[0])

D:\Users\sampl\Anaconda3\lib\site-packages\sympy\core\sympify.py in _sympify(a)
    415 
    416     """
--> 417     return sympify(a, strict=True)
    418 
    419 

D:\Users\sampl\Anaconda3\lib\site-packages\sympy\core\sympify.py in sympify(a, locals, convert_xor, strict, rational, evaluate)
    337 
    338     if strict:
--> 339         raise SympifyError(a)
    340 
    341     if iterable(a):

SympifyError: SympifyError: array([[1353.5478432 - 4.955328*Tinfhout,
        1363.55860683636 - 4.955328*Tinfhout,
        1373.56937047273 - 4.955328*Tinfhout, ...,
        2324.59191592727 - 4.955328*Tinfhout,
        2334.60267956364 - 4.955328*Tinfhout,

Here is my code:

from pylab import *
from random import *
from mpl_toolkits import mplot3d
import pandas as pd
from scipy.optimize import fsolve
import sympy as sy

mdoth = 0.004916
cph = 1008
nsh = .598
hh= 86.68
Ash = .02
n=127
alpha = .00041427
rho = .002129
k=3.041
Le = .0025
Ae = .000001
re = rho * Le/Ae
Ke = k * Ae/Le
nsc = .674
hc = 87.68
Asc = .016
rL = re
mdotc = .004542
cpc = 1007




II, Qc, Qh, Tc, Th, Tinfcout, Tinfhout = symbols('II, Qc, Qh, Tc, Th, Tinfcout, Tinfhout')

def solver(TH, TC):
    Tinfhin = TH +273.15
    Tinfcin = TC + 273.15
    sols = sy.nsolve(  (Eq(Qh,mdoth * cph * (Tinfhin - Tinfhout) ),
          Eq(Qh,nsh * hh * Ash * ((Tinfhin + Tinfhout)/2 - Th)),
          Eq(Qh,n * (alpha * II * Th - 0.5 * (II**2) * ree + (Ke * (Th-Tc)))), 
          Eq(Qc,n * (alpha * II * Tc + 0.5 * (II**2) * ree + (Ke * (Th-Tc)))),
          Eq(Qc,nsc * hc * Asc * (Tc - (Tinfcin + Tinfcout)/2) ),
          Eq(Qc,mdotc * cpc * (Tinfcout - Tinfcin) ),
          Eq(II,(alpha * (Th - Tc))/(rL + ree) )),
      (II, Qc, Qh, Tc, Th, Tinfcout, Tinfhout), (1,5,5,300,300,330,330) )
    result = sols[0]
    return(result)     


hin = linspace(0,200,100)
cin = linspace(0,200,100)
X, Y = meshgrid(hin,cin)
Z = solver(X,Y)

ax = axes(projection='3d')
ax.set_xlabel("TC")
ax.set_ylabel("Ambient")
ax.set_zlabel("Voltage")
ax.plot_surface(X, Y, Z, rstride=1, cstride=1, cmap = 'plasma')
ax.view_init(0, 180)'''

What is the best solution to this problem?

Answer 1

Do not use these, it is a bad practice. And please use ìmport matplotlib.pyplot.

from pylab import *
from random import *

The improved code is:

import matplotlib
import matplotlib.pyplot as plt
import random
# from mpl_toolkits import mplot3d
# import pandas as pd
from scipy.optimize import fsolve
import sympy as sy
import numpy as np

mdoth = 0.004916
cph = 1008
nsh = .598
hh= 86.68
Ash = .02
n=127
alpha = .00041427
rho = .002129
k=3.041
Le = .0025
Ae = .000001
ree = rho * Le/Ae
Ke = k * Ae/Le
nsc = .674
hc = 87.68
Asc = .016
rL = ree
mdotc = .004542
cpc = 1007




II, Qc, Qh, Tc, Th, Tinfcout, Tinfhout = sy.symbols('II, Qc, Qh, Tc, Th, Tinfcout, Tinfhout')

def solver(_TH, _TC):
    rv = []
    for TH,TC in zip(_TH, _TC):
        TH = TH[0]
        TC = TC[0]
        Tinfhin = TH +273.15
        Tinfcin = TC + 273.15
        sols = sy.nsolve((sy.Eq(Qh,mdoth * cph * (Tinfhin - Tinfhout) ),
            sy.Eq(Qh,nsh * hh * Ash * ((Tinfhin + Tinfhout)/2 - Th)),
            sy.Eq(Qh,n * (alpha * II * Th - 0.5 * (II**2) * ree + (Ke * (Th-Tc)))), 
            sy.Eq(Qc,n * (alpha * II * Tc + 0.5 * (II**2) * ree + (Ke * (Th-Tc)))),
            sy.Eq(Qc,nsc * hc * Asc * (Tc - (Tinfcin + Tinfcout)/2) ),
            sy.Eq(Qc,mdotc * cpc * (Tinfcout - Tinfcin) ),
            sy.Eq(II,(alpha * (Th - Tc))/(rL + ree) )),
        (II, Qc, Qh, Tc, Th, Tinfcout, Tinfhout), (1,5,5,300,300,330,330) )
        rv.append(sols[0])
    return rv


hin = np.linspace(0, 200, 20)
cin = np.linspace(0, 200, 20)
X, Y = np.meshgrid(hin,cin)
Z = solver(X,Y)

ZZ = []

for _ in range(0, len(Z)):
    ZZ.append(Z)
    
ZZ = np.array(ZZ, dtype='float')

fig, ax = plt.subplots(figsize=(8, 8), subplot_kw={"projection": "3d"})
ax.plot_surface(X, Y, ZZ, rstride=1, cstride=1, cmap = 'plasma', antialiased=False)

ax.set_xlabel("TC")
ax.set_ylabel("Ambient")
ax.set_zlabel("Voltage")
ax.view_init(0, 180)
fig.tight_layout()
plt.show()

The figure is

I have used the function made by @smichr.

Answer 2

One must always need more care when using multiple packages since idioms in one are not always applicable in the other. SymPy is telling you that it doesn't know what to do with the array object. I am thinking you will need to unpack the arrays elements one at a time to solve and build up the solution vector. *And also change the variable name re to ree everywhere:

def solver(_TH, _TC):
    rv = []
    for TH,TC in zip(_TH, _TC):
        TH = TH[0]
        TC = TC[0]
        Tinfhin = TH +273.15
        Tinfcin = TC + 273.15
        sols = sy.nsolve(  (Eq(Qh,mdoth * cph * (Tinfhin - Tinfhout) ),
            Eq(Qh,nsh * hh * Ash * ((Tinfhin + Tinfhout)/2 - Th)),
            Eq(Qh,n * (alpha * II * Th - 0.5 * (II**2) * ree + (Ke * (Th-Tc)))), 
            Eq(Qc,n * (alpha * II * Tc + 0.5 * (II**2) * ree + (Ke * (Th-Tc)))),
            Eq(Qc,nsc * hc * Asc * (Tc - (Tinfcin + Tinfcout)/2) ),
            Eq(Qc,mdotc * cpc * (Tinfcout - Tinfcin) ),
            Eq(II,(alpha * (Th - Tc))/(rL + ree) )),
        (II, Qc, Qh, Tc, Th, Tinfcout, Tinfhout), (1,5,5,300,300,330,330) )
        rv.append(sols[0])
    return(rv)