fsolve gives weird answers

Question

I want to use fsolve to numerically find roots of a nonlinear transcendent equation.

The following code does this job.

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

kappa = 0.1
tau = 90

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

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()

I can double check the solutions graphically by plotting the two graphs f1(x)=x and f2(x)=k * sin(-x * tau), which i also included in the code.
fsolve gives me some wrong answers, without throwing any errors or convergence problems.

The Problem is, that I would like to automatize the procedure for varying kappa and tau, without me checking which answers are wrong and which are right. But with wrong answers as output, i can't use this method. Is there any other method or an option I can use, to be on the safe side?

Thanks for the help.

RenZhen95 RenZhen95 · Accepted Answer · 2020-04-24T23:37:17

I couldn't run your code, there were errors and anyway, according to the documentation on scipy.fsolve, you're supposed to add an initial guess as the second input argument, not a range as what you did there fsolve(equation, x0, (kappa,tau))

You could of course however pass this in a loop, looping for every value in the array np.linspace(0.5, 0.5, 25). Although I do not understand what you are trying to achieve by varying kappa and tau, but if I take it that for those given parameters, you are interested in looking for the roots, here's how I would do it.

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]

Probably not the smartest way to do it (assuming I even understood you correctly), but perhaps you have an idea now.

fsolve gives weird answers

1 Answers