Решить нелинейное уравнение в Python

Я пытаюсь найти фундаментальную ТЕ-моду диэлектрического волновода. Я пытаюсь решить эту проблему, чтобы вычислить две функции и попытаться найти их пересечение на графике. Однако у меня возникли проблемы с получением точки пересечения на графике. Мой код:

def LHS(w):
    theta = 2*np.pi*1.455*10*10**(-6)*np.cos(np.deg2rad(w))/(900*10**(-9))
    if(theta>(np.pi/2) or theta < 0):
        y1 = 0
    else:
        y1 = np.tan(theta)
    return y1

def RHS(w):
    y = ((np.sin(np.deg2rad(w)))**2-(1.440/1.455)**2)**0.5/np.cos(np.deg2rad(w))
    return y

x = np.linspace(80, 90, 500)

LHS_vals = [LHS(v) for v in x]
RHS_vals = [RHS(v) for v in x]


# plot
fig, ax = plt.subplots(1, 1, figsize=(6,3))
ax.plot(x,LHS_vals)
ax.plot(x,RHS_vals)
ax.legend(['LHS','RHS'],loc='center left', bbox_to_anchor=(1, 0.5))

plt.ylim(0,20)
plt.xlabel('degree')
plt.ylabel('magnitude')
plt.show()

У меня получился такой сюжет: graph

Точка пересечения составляет около 89 градусов, однако мне сложно вычислить точное значение x. Я пробовал fsolve, решить, чтобы найти решение, но все равно напрасно. Кажется, невозможно распечатать решение, если это не единственное решение. Можно ли найти решение только в том случае, если x находится в определенном диапазоне? Может ли кто-нибудь дать мне здесь какое-нибудь предложение? Спасибо!

изменить: уравнение выглядит так (m = 0): введите описание изображения здесь

и я пытаюсь решить тета здесь, найдя точку пересечения

edit: Один из способов, который я пробовал, таков:

from scipy.optimize import fsolve
def f(wy):
   w, y = wy
   z = np.array([y - LHS(w), y - RHS(w)])
   return z

fsolve(f,[85, 90])

Однако это дает мне неправильный ответ. Я тоже пробовал что-то вроде этого:

import matplotlib.pyplot as plt

x = np.arange(85, 90, 0.1)
f = LHS(x)
g = RHS(x)

plt.plot(x, f, '-')
plt.plot(x, g, '-')

idx = np.argwhere(np.diff(np.sign(f - g)) != 0).reshape(-1) + 0
plt.plot(x[idx], f[idx], 'ro')
print(x[idx])
plt.show()

Но он показывает: ValueError: истинное значение массива с более чем одним элементом неоднозначно. Используйте a.any () или a.all ()


python numpy scipy equation-solving
person tsen0406    schedule 12.02.2018    source источник
comment
Можете ли вы опубликовать уравнение, которое пытаетесь решить? Кроме того, как вы ищете перекресток?   -  person Alex Dubrovsky    schedule 12.02.2018
comment
Обязательно пометьте вопрос в соответствии с актуальной проблемой и покажите код актуальной проблемы. В вопросе нет попытки использовать fsolve или что-то подобное. (С другой стороны, проблема совершенно не связана с matplotlib.)   -  person ImportanceOfBeingErnest    schedule 12.02.2018
comment
Во-первых, ваша LHS является прерывистой функцией, fsolve предполагает плавный аргумент. Во-вторых, ваш RHS отличается от того, что вы показываете на своем участке. В-третьих, ваша RHS не является реальной функцией, т.е. аргумент sqrt(x) становится отрицательным, когда w приближается к значению 82.   -  person Dmitri Chubarov    schedule 12.02.2018
comment
Общее замечание из моего опыта работы с нелинейными уравнениями - тщательно исследуйте уравнение, прежде чем решать его численно, особенно для функций с разрывами и таких вещей, как отрицательные аргументы sqrt(x).   -  person Alex Dubrovsky    schedule 12.02.2018

Ответы (2)


arrow_upward
2
arrow_downward

Что-то быстрое и (очень) грязное, которое, похоже, работает (по крайней мере, оно дало значение тета ~ 89 для ваших параметров) - добавьте в свой код перед цифрой после строки RHS_vals = [RHS(v) for v in x] следующее:

# build a list of differences between the values of RHS and LHS
diff = [abs(r_val- l_val) for r_val, l_val in zip(RHS_vals, LHS_vals)]

# find the minimum of absolute values of the differences
# find the index of this minimum difference, then find at which angle it occured
min_diff = min(diff)
print "Minimum difference {}".format(min_diff)
print "Theta = {}".format(x[diff.index(min_diff)])

Я смотрел в диапазоне 85-90:

x = np.linspace(85, 90, 500)
person Alex Dubrovsky    schedule 12.02.2018

arrow_upward
3
arrow_downward

Во-первых, вам нужно убедиться, что функция действительно может обрабатывать массивы numpy. Несколько вариантов определения кусочных функций показаны в Построение дискретного распределения с использованием np.linspace () . Например.

def LHS(w):
    theta = 2*np.pi*1.455*10e-6*np.cos(np.deg2rad(w))/(900e-9)
    y1 = np.select([theta < 0, theta <= np.pi/2, theta>np.pi/2], [np.nan, np.tan(theta), np.nan])
    return y1

Это уже позволяет использовать второй подход, нанося точку с индексом, наиболее близким к минимуму разницы между двумя кривыми.

import numpy as np
import matplotlib.pyplot as plt

def LHS(w):
    theta = 2*np.pi*1.455*10e-6*np.cos(np.deg2rad(w))/(900e-9)
    y1 = np.select([theta < 0, theta <= np.pi/2, theta>np.pi/2], [np.nan, np.tan(theta), np.nan])
    return y1

def RHS(w):
    y = ((np.sin(np.deg2rad(w)))**2-(1.440/1.455)**2)**0.5/np.cos(np.deg2rad(w))
    return y

x = np.linspace(82.1, 89.8, 5000)
f = LHS(x)
g = RHS(x)

plt.plot(x, f, '-')
plt.plot(x, g, '-')

idx = np.argwhere(np.diff(np.sign(f - g)) != 0).reshape(-1) + 0
plt.plot(x[idx], f[idx], 'ro')
plt.ylim(0,40)
plt.show()

введите здесь описание изображения

Затем можно также использовать scipy.optimize.fsolve, чтобы получить фактическое решение.

idx = np.argwhere(np.diff(np.sign(f - g)) != 0)[-1]

from scipy.optimize import fsolve

h = lambda x: LHS(x)-RHS(x)
sol = fsolve(h,x[idx])

plt.plot(sol, LHS(sol), 'ro')
plt.ylim(0,40)
plt.show()
person ImportanceOfBeingErnest    schedule 13.02.2018