Почему корни этой функции, найденные scipy.optimize с использованием метода деления пополам, не равны нулю? - PullRequest
Новая
       105

Почему корни этой функции, найденные scipy.optimize с использованием метода деления пополам, не равны нулю?

0 голосов
/ 09 июля 2020

Я пытаюсь решить уравнение, которое можно решить только численными методами. Проблема, которую я обнаружил, заключается в том, что при использовании метода пополам для нахождения root этого уравнения функция не оценивает ноль в некоторых найденных корнях.

В частности, я пытаюсь чтобы найти корни func, определенного как: 

import numpy as np
from scipy.optimize import bisect
import matplotlib.pyplot as plt

d = 0.01

c_L = 6420
c_T = 3040

def func(y, x):
    w = 2*np.pi*(x/d)
    k = w/y

    p = np.sqrt((w/c_L)**2 - k**2, dtype = np.complex128)
    q = np.sqrt((w/c_T)**2 - k**2, dtype = np.complex128)

    func = np.tan(q*d/2)/q + (4*(k**2)*p*np.tan(p*d/2))/(q**2 - k**2)**2

    return np.real(func)

Функция является непрерывной, поэтому я могу применить метод деления пополам. Предположим, я фиксирую значение x = 1600 и решаю для y. В интервале [4900, 5000] функция оценивается с противоположными знаками: 

p = func(4900, 1600)
q = func(5000, 1600)

print(np.sign(p))
print(np.sign(q))

Результат: -1 и 1. Применение метода деления пополам к этому интервалу и вычисление root, найденного в функции, дает: 

bisection_pq = bisect(func, 4900, 5000, args = (1600,))

print(func(bisection_pq, 1600))

Результат: 9.194034422677078e-17, что близко к нулю, так что это допустимое решение. Теперь функция снова меняет знак с интервалом [9700, 9800]. Выполнение той же процедуры, что и раньше: 

s = func(9700, 1600)
t = func(9800, 1600)

print(np.sign(s))
print(np.sign(t))

bisection_st = bisect(func, 9700, 9800, args = (1600,))

print(func(bisection_st, 1600))

Теперь на выходе будет -8314070119023.9795, что не близко к нулю, поэтому это не root. Я действительно не понимаю, что происходит.

В конечном итоге я пытаюсь построить график решения уравнения как функции от x. Если я исправлю кратные значения x и решу для y так же, как и раньше, сценарий будет выглядеть так: 

nmodes = 6

x_arr = np.arange(0, 10000, 100)
y_arr = np.array([[start, start + 100] for start in range(0, 15000, 100)])

result = np.zeros((len(x_arr), nmodes))
result_evaluated_at_bisect = np.zeros((len(x_arr), nmodes))

for i, x in enumerate(x_arr):

    print(str(i) + '/100', x, sep = ' - ')

    j = 0

    for y in y_arr:

        a = func(y[0], x)
        b = func(y[1], x)

        if j < nmodes:
            if not np.isnan(a) and not np.isnan(b):
                if np.sign(a) != np.sign(b):
                    bisection = bisect(func, y[0], y[1], args = (x,))
                    result[i][j] = bisection
                    result_evaluated_at_bisect[i][j] = func(bisection, x)
                    j += 1

plt.scatter(np.tile(x_arr, (3,1)).T, result[:, 1::2], label = 'Incorrect')
plt.scatter(np.tile(x_arr, (3,1)).T, result[:, ::2], label = 'Correct')

plt.xlabel('x')
plt.ylabel('y')

plt.ylim([200, 15000])
plt.xlim([0, 10000])

plt.legend()

plt.show()

Результат этого сценария будет выглядеть как this . Точки, помеченные как «Неверно», не равны нулю при оценке функции в root, найденном с помощью optimize.bisect. Те, которые помечены как «Правильные», должны быть единственными ожидаемыми результатами (я проверил эти значения с помощью проверенного программного обеспечения). Что-то не так в моей процедуре?

Заранее спасибо.

...