I am quite new to python and I have been trying to plot this in a few different ways. If I try using np.vectorize, it crashes.

so i wrote this code, which is giving me the error in the title:

import math
import numpy as np
from scipy import optimize
from scipy import interpolate
from scipy import integrate
import matplotlib.pyplot as plt
from sympy import Symbol
from sympy.solvers import solve
import sympyfig = plt.figure(2, figsize=(10, 7))
ax = fig.gca(projection='3d')
x = np.linspace(-1, 1, num=20)
y = np.linspace(-1, 1, num=20)
X, Y = np.meshgrid(x, y)#we define a function where we use sympy.solvers.solve and return both values of z
def surfacez(x, y):z = Symbol('z')surface = solve(x + 2 * y + z + math.e ** (2 * z), z)return surfaceZ = surfacez(X,Y)#we plot for both values of z
surf = ax.plot_wireframe(X, Y, Z, linewidth=1, antialiased=False, color='b')
plt.xlabel('x')
plt.ylabel('y')
ax.set_zlabel('z')
plt.title('Quadratic surface')

with full traceback:

Traceback (most recent call last):

  File "C:\Users\natty\Desktop\LU\MATB21\crash.py", line 39, in <module>Z = surfacez(X,Y)File "C:\Users\natty\Desktop\LU\MATB21\crash.py", line 32, in surfacezsurface = solve(x + 2 * y + z + math.e ** (2 * z), z)File "C:\Users\natty\AppData\Local\Programs\Spyder\pkgs\sympy\solvers\solvers.py", line 1096, in solvesolution = _solve_system(f, symbols, **flags)File "C:\Users\natty\AppData\Local\Programs\Spyder\pkgs\sympy\solvers\solvers.py", line 1730, in _solve_systemi, d = _invert(g, *symbols)File "C:\Users\natty\AppData\Local\Programs\Spyder\pkgs\sympy\solvers\solvers.py", line 3109, in _invertindep, dep = lhs.as_independent(*symbols)AttributeError: 'ImmutableDenseNDimArray' object has no attribute 'as_independent'

i do not know what this error means, could you suggest a way to fix or the origin of this error?

In [66]: z
Out[66]: z
In [67]: X,Y=np.meshgrid([1,2,3],[1,2])
In [68]: X
Out[68]: 
array([[1, 2, 3],[1, 2, 3]])
In [69]: Y
Out[69]: 
array([[1, 1, 1],[2, 2, 2]])
In [71]: expr = X +2 * Y + z + E  ** (2*z)
In [72]: expr
Out[72]: 
array([[z + exp(2*z) + 3, z + exp(2*z) + 4, z + exp(2*z) + 5],[z + exp(2*z) + 5, z + exp(2*z) + 6, z + exp(2*z) + 7]],dtype=object)

Looks like solve has converted this numpy array to sympy Array:

In [74]: Array(expr)
Out[74]: 
⎡     2⋅z           2⋅z           2⋅z    ⎤
⎢z + ℯ    + 3  z + ℯ    + 4  z + ℯ    + 5⎥
⎢                                        ⎥
⎢     2⋅z           2⋅z           2⋅z    ⎥
⎣z + ℯ    + 5  z + ℯ    + 6  z + ℯ    + 7⎦
In [75]: type(_)
Out[75]: sympy.tensor.array.dense_ndim_array.ImmutableDenseNDimArray

I haven't explored solve that much, but evidently this isn't one of the kinds of systems of expressions that it can handle.

For scalar values

In [80]: X,Y=1,2
In [83]: expr = X +2 * Y + z + E  ** (2*z)
In [84]: expr
Out[84]: 2⋅z    
z + ℯ    + 5
In [85]: solve(expr,z)
Out[85]: 
⎡      ⎛   -10⎞⎤
⎢     W⎝2⋅ℯ   ⎠⎥
⎢-5 - ─────────⎥
⎣         2    ⎦

and with symbols x and y:

In [86]: expr = x +2 * y + z + E  ** (2*z)
In [87]: expr
Out[87]: 2⋅z
x + 2⋅y + z + ℯ   
In [88]: solve(expr,z)
Out[88]: 
⎡            ⎛   -2⋅x - 4⋅y⎞⎤
⎢           W⎝2⋅ℯ          ⎠⎥
⎢-x - 2⋅y - ────────────────⎥
⎣                  2        ⎦
In [89]: sol=_[0]
In [90]: sol.subs({x:1, y:2})
Out[90]: ⎛   -10⎞W⎝2⋅ℯ   ⎠
-5 - ─────────2    

We could try creating a python function from that:

In [91]: f = lambdify((x,y),sol)
In [92]: help(f)Help on function _lambdifygenerated:_lambdifygenerated(x, y)Created with lambdify. Signature:func(x, y)Expression:-x - 2*y - LambertW(2*exp(-2*x - 4*y))/2Source code:def _lambdifygenerated(x, y):return (-x - 2*y - 1/2*lambertw(2*exp(-2*x - 4*y)))

I'm not sure about that LambertW. Obviously there's nothing in the stock python/numpy like that. scipy maybe.

OK, scipy does have it: https://docs.scipy.org/doc/scipy/reference/generated/scipy.special.lambertw.html

In [95]: from scipy.special import lambertw as lambertW
In [96]: f(1,2)
Out[96]: (-5.000045395808017+0j)

Compare that with sympy's own numeric solution:

In [97]: sol.subs({x:1, y:2})
Out[97]: ⎛   -10⎞W⎝2⋅ℯ   ⎠
-5 - ─────────2    
In [98]: N(sol.subs({x:1, y:2}))
Out[98]: -5.00004539580802

and with array arguments, the lambdified function:

In [99]: X,Y=np.meshgrid([1,2,3],[1,2])
In [100]: f(X,Y)
Out[100]: 
array([[-3.00246655+0.j, -4.00033524+0.j, -5.0000454 +0.j],[-5.0000454 +0.j, -6.00000614+0.j, -7.00000083+0.j]])

sympy.lambdify is one of the best tools for using sympy expressions in numpy functions. It can't handle all cases, since the translation is largely lexical, without deep understanding of numpy. There have been SO about such cases.