physical problem

[Sofiaphys]

I am currently writing a code that solves a system of 4 ODES and then uses the solutions to calculate some quantities (W1, W2) in the code, for different values of two parameters n, sigma. The goal is to find the values of the parameters n, sigma for wich the quantites W1, W2 are equal to zero and plot them in a plot n=n(sigma). Even though the code runs without errors, it has a warning about the ODE system that i do not understand and also it runs very slowly for a big number of values for the n and sigma's. What's the problem with it? Also the last time i ran it, it did not plot anything.

import numpy as np
import matplotlib.pyplot as plt
import scipy as sc
from scipy.integrate import odeint
from scipy.interpolate import CloughTocher2DInterpolator

m=[]
X=[]
Y=[]
Z=[]
Z2=[]
M=1

for sigma in np.arange(-5, 5, 0.1):
    for n in np.arange(-5, 5, 0.1):
        
            def dSdx(S, x):
                y1, y2, p1, p2 = S
                s=sigma-np.tanh(x)
                o=(s**2+n**2)**2
                return[(M**2)*(((s**2-n**2)/o)-1)*p1+2*(M**2)*s*n*p2/o, (M**2)*(((s**2-n**2)/o)-1)*p2-2*(M**2)*s*n*p1/o, (s**2-n**2)*y1/(M**2)-2*s*n*y2/(M**2), (s**2-n**2)*y2/(M**2)+2*s*n*y1/(M**2)]

            y1_0=0
            y2_0=0 
            p1_0=1
            p2_0=0
        
            S_0=(y1_0, y2_0, p1_0, p2_0)
            x=np.linspace(-1, 1, 100)
            sol=odeint(dSdx, S_0, x)
        
            y1_1=sol.T[0][-1]
            y2_1=sol.T[1][-1]
            p1_1=sol.T[2][-1]
            p2_1=sol.T[3][-1]
        
            W2=y1_1*p2_1-y2_1*p1_1
            W1=y1_1*p1_1+y2_1*p2_1
            
            k=[n, sigma, W1, W2]
           
            m.append(k)
            X.append(n)
            Y.append(sigma)
            Z.append(W1)
            Z2.append(W2)


interp_func=CloughTocher2DInterpolator(list(zip(X, Y)), Z)
interp_func2=CloughTocher2DInterpolator(list(zip(X, Y)), Z2)
x_values=np.arange(-5, 5, 0.05)
y_values=np.arange(-5, 5, 0.05)
x_zero_contour=[]
y_zero_contour=[]
x_zero_contour2=[]    
y_zero_contour2=[]

for x_val in x_values:
    for y_val in y_values:
        if np.isclose(interp_func(x_val, y_val), 0, atol=1e-4):
            x_zero_contour.append(x_val)
            y_zero_contour.append(y_val)
            
for x_val in x_values:
    for y_val in y_values:
        if np.isclose(interp_func2(x_val, y_val), 0, atol=1e-4):
            x_zero_contour2.append(x_val)
            y_zero_contour2.append(y_val)


plt.plot(y_zero_contour, x_zero_contour, label='solution path')
plt.plot(y_zero_contour2, x_zero_contour2, label='conjugate path')
plt.xlabel('sigma')
plt.ylabel('n')
plt.legend()

buran write Nov-22-2023, 08:08 PM:
Please, use proper tags when post code, traceback, output, etc. This time I have added tags for you.
See BBcode help for more info.

[Gribouillis]

I added at line 16

print(sigma, n)

and I ran the code with python -W error. The output shows where the warning occurs

Output:
...
-0.7000000000000153 -0.6000000000000156
-0.7000000000000153 -0.500000000000016
-0.7000000000000153 -0.40000000000001634
-0.7000000000000153 -0.3000000000000167
-0.7000000000000153 -0.20000000000001705
-0.7000000000000153 -0.10000000000001741
-0.7000000000000153 -1.7763568394002505e-14
Traceback (most recent call last):
  File "/home/eric/Projets/Scratch/2023-01/paillasse/pf/sophiaphys.py", line 31, in <module>
    sol=odeint(dSdx, S_0, x)
  File "/usr/lib/python3/dist-packages/scipy/integrate/_odepack_py.py", line 247, in odeint
    warnings.warn(warning_msg, ODEintWarning)
scipy.integrate._odepack_py.ODEintWarning: Excess work done on this call (perhaps wrong Dfun type). Run with full_output = 1 to get quantitative information.
λ 

So the error occurs at some point where n is almost equal to 0. This may help find why the solver fails at this point.