fsolve дает странные ответы

Question

Я хочу использовать fsolve для численного поиска корней нелинейного трансцендентного уравнения.

Следующий код выполняет эту работу.

import numpy as np
from scipy.optimize import fsolve
import matplotlib.pyplot as plt</p>

<p>kappa = 0.1
tau = 90</p>

<p>def equation(x, * parameters):
    kappa,tau = parameters
    return -x + kappa * np.sin(-tau*x)</p>

<code>x = np.linspace(-0.5,0.5, 35)
roots = fsolve(equation,x, (kappa,tau))


x_2 = np.linspace(-1.5,1.5,1500)
plt.plot(x_2 ,x_2 )
plt.plot(x_2 , kappa*np.sin(-x_2 *tau))
plt.scatter(x, roots)
plt.show()
</code>

Я могу дважды проверить решения графически, построив два графика f1 (x) = x и f2 (x) = k * sin (-x * tau), которые я также включил в код. fsolve дает мне несколько неправильных ответов, без каких-либо ошибок или проблем сходимости.

Проблема в том, что я хотел бы автоматизировать процедуру изменения каппа и тау, не проверяя, какие ответы неправильные, а какие - правильные. , Но с неправильными ответами в качестве вывода, я не могу использовать этот метод. Есть ли какой-либо другой способ или вариант, который я могу использовать, чтобы быть в безопасности?

Спасибо за помощь.

Answer 1

Я не смог запустить ваш код, были ошибки и в любом случае, согласно документации на scipy.fsolve , вы должны добавить начальное предположение в качестве второго входного аргумента, а не диапазона как то, что вы сделали там fsolve(equation, x0, (kappa,tau))

Однако вы, конечно, можете передать это в al oop, выполняя цикл для каждого значения в массиве np.linspace(0.5, 0.5, 25). Хотя я не понимаю, чего вы пытаетесь достичь, изменяя каппу и тау, но если я возьму это для этих заданных параметров, вы заинтересованы в поиске корней, вот как я бы это сделал.

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

# Take it as it is
kappa   = 0.1
tau     = 90

def equation(x, parameters):
    kappa,tau = parameters
    return -x + kappa * np.sin(-tau*x)

# Initial guess of x = -0.1
SolutionStack   = []
x0              = -kappa
y               = fsolve(equation, x0, [kappa, tau])
SolutionStack.append(y[0])

y               = fsolve(equation, SolutionStack[-1], [kappa, tau])
SolutionStack.append(y[0])
deltaY = SolutionStack[-1] - SolutionStack[0]

# Define tolerance
tol = 5e-4
while ((SolutionStack[-1] <= kappa) and (deltaY <= tol)):
    y = fsolve(equation, SolutionStack[-1], [kappa, tau])
    SolutionStack.append(y[0])
    deltaY = SolutionStack[-1] - SolutionStack[-2]
    # Obviously a little guesswork is involved here, as it pertains to 0.07
    if deltaY <= tol:
        SolutionStack[-1] = SolutionStack[-1] + 0.07

# Obtaining the roots
Roots = []
Roots.append(SolutionStack[0])
for i in range(len(SolutionStack)-1):
    if (SolutionStack[i+1] - SolutionStack[i]) <= tol:
        continue
    else:
        Roots.append(SolutionStack[i+1]

Возможно, это не самый умный способ сделать это (если я вас правильно понял), но, возможно, у вас есть идея сейчас.