Во-первых, numpy и sympy не очень хорошо сочетаются. Sympy не знает о функциях numpy, а numpy не понимает символов sympy. Итак, лучше всего выполнить все вычисления Symboli c в Sympy. И, наконец, lambdify может преобразовать выражения sympy в numpy функции.
Вы можете немного помочь sympy, указав больше информации о переменных, например, зная, что positive=True может помочь упростить квадратные корни .
Что касается черчения, то есть некоторые дополнения к sympy для быстрого показа сюжета. Но они довольно ограничены в возможностях, особенно если вы хотите 3D-графики. Для трехмерных графиков matplotlib - это пакет goto, а numpy - предпочтительный способ представления функций и данных. Тем более, что у ваших функций огромный диапазон чисел.
Обратите внимание, что solve(est) без дополнительных уточнений решает для d как функцию rep. Назовите его как solve(est, rep), если вам нужно rep как функция d.
В приведенном ниже коде все функции numpy во время пути symboli c изменены на функции sympy. Затем lambdify используется для создания трехмерного графика поверхности с помощью matplotlib. Обратите внимание, что lambdify не может обрабатывать уравнение напрямую, но его нужно переписать как left-hand-side - right-hand-side.
(в этом коде используется Python 3. Единственное отличие состоит в том, что для печати нужны скобки. Если нет важных причина использования Python 2, так как с 2020 года очень рекомендуется перейти на Python 3.)
import sympy as sym
nc = 1.772
c = 3e10
L = 0.3
R = 3.0
lbd = 775e-7
d = sym.Symbol('d', real=True, positive=True, nonzero=True)
rep = sym.Symbol('rep', real=True, positive=True, nonzero=True)
M1 = sym.Matrix([[1, (L / 2)], [0, 1]]) # espessura_cristal
M2 = sym.Matrix([[1, 0], [0, (1 / nc)]]) # ar_cristal
M3 = sym.Matrix([[1, (d / 2) - L / 2], [0, 1]]) # dis_cristal_esp
M4 = sym.Matrix([[1, 0], [(-2 / R), 1]]) # espelho
M5 = sym.Matrix([[1, (c / rep) - d], [0, 1]]) # comp_cavidade
M6 = sym.Matrix([[1, 0], [0, nc]]) # cristal_ar
M = sym.Matrix(M1 * M2 * M3 * M4 * M5 * M4 * M3 * M6 * M1)
est = sym.Eq(((M[0, 0] + M[1, 1]) ** 2), 4.0 )
sol = sym.solve(est, rep)
print("sol:", sol)
print("est:", est)
# sym.plot_implicit(est, (rep, 0.1, 1), (d, 1, 10))
est_np = sym.lambdify((rep, d), (est.lhs - est.rhs).simplify(), "numpy")
import numpy as np
from matplotlib import pyplot as plt
from mpl_toolkits import mplot3d
rep_np, d_np = np.meshgrid(np.linspace(0.1, 1, 200), np.linspace(1, 10, 200))
fig = plt.figure()
ax = plt.axes(projection='3d')
Z = est_np(rep_np, d_np)
ax.plot_surface(rep_np, d_np, Z, rstride=1, cstride=1, cmap='plasma', edgecolor='none')
ax.set_xlabel('rep')
ax.set_ylabel('d')
ax.set_zlabel('est')
plt.show()
Во время пребывания в Однако вы можете построить 4 решения следующим образом (после вызова sol = sym.solve(est, d)):
plots = []
for i, s in enumerate(sol):
plots.append(sym.plot(s, xlim=(.001, 1), xlabel='rep', ylabel='d', title='solution '+ str(i+1), show=False))
sym.plotting.PlotGrid(2, len(sol)//2, *plots)
Первые два решения выглядят как очень похоже на sym.plot_implicit вопроса. Последние два решения получают огромные значения для d. Может быть, это не реальные решения, но из-за проблем с округлением?