Problems to solve non-linear coupled dgl in python (numerical)

Question

I have the following problem. I would like to solve the following DGL system numerically in python:

$$ \begin{align*} &g'' + \frac{2}{r} g' - \frac{2}{r^2} g - \frac{3}{r} e g^2 - e^2 g^3 - \frac{e}{2r} f^2 - \frac{e^2}{4} f^2g = 0 \ &f'' + \frac{2}{r} f' - \frac{2}{r^2} f - \frac{2e}{r} fg - \frac{e^2}{2} fg^2 + \nu^2 f - \lambda f^3 = 0 \end{align*} $$

The procedure should be the same as always. I have 2 coupled differential equations of the 2nd order and I use the substitution g' = v and f' = u to create four coupled differential equations of the 1st order.

\begin{align*} \frac{dg}{dr} &= v \ \frac{dv}{dr} &= -\frac{2}{r}v + \frac{2}{r^2}g - \frac{3}{r}g^2 - g^3 - \frac{2}{r}f^2 - fg \ \frac{df}{dr} &= u \ \frac{du}{dr} &= -\frac{2}{r}u - \frac{2}{r^2}f - \frac{2}{r}gf - \frac{1}{2}gf^2 + \frac{1}{r^2}f - \frac{\alpha}{r^2}f^3 \end{align*}

Now I should be able to calculate them numerically or not? I have used the packet scipy and tried various solution methods (BDF, RK45, LSODA) but I am not getting a solution, only 0. I have also written my own RK4, but I have the same problem here too. I have also tried to set boundary conditions in the infinite (f(inf) = 1, g(inf) = 0 ) nothing works.

Here is my code:

from scipy.integrate import solve_ivp
import numpy as np
import matplotlib.pyplot as plt

# Constants
e = 1
mu = 1
lambd = 1

# Initial second-order system

# g'' + 2/r g'-2/r^2 g - 3/r eg^2- e^2 g^3 - e/(2r) f^2 - e^2/4 f^2g = 0
# f'' + 2/r f' - 2/r^2 f - 2e/r fg - e^2/2 fg^2 + mu^2 f - lambda f^3 = 0

# Substitution
# g' = v
# v' = g''
# f' = u
# u' = f''

# Equivalent set of first-order systems
def dSdr(r, S):
    g, v, f, u = S
    dgdr = v
    dvdr = -2 / r * v - 2 / r ** 2 * g + 3 / r * e * g ** 2 + e ** 2 * g ** 3 + e / (
                2 * r) * f ** 2 + e ** 2 / 4 * f ** 2 * g
    dfdr = u
    dudr = -2 / r * u + 2 / r ** 2 * f - 2 * e / r * f * g + e ** 2 / 2 * f * g ** 2 - mu ** 2 * f + lambd * f ** 3

    # Correction of boundary conditions for large r
    if r > 10000:
        dgdr = 0  # g(infinity) = 0
        dfdr = 1  # f(infinity) = 1

    return [dgdr, dvdr, dfdr, dudr]

# Equivalent set of first-order systems (Scaled)
def dSdr_scaled(r, S):
    g, v, f, u = S
    alpha = 1
    dgdr = v
    dvdr = -2/r * v + 2/r**2 * g - 3 / r * g**2 - g**3 - 2/r * f**2 - f*g
    dfdr = u
    dudr = -2 / r * u * 2/r**2 * f - 2/r * g * f - 1/2 * g * f**2 + 1/r**2 * f - alpha/r**2 * f**3
    return [dgdr, dvdr, dfdr, dudr]

# Initial conditions
g0 = 0
v0 = 0
f0 = 0
u0 = 0

S0 = [g0, v0, f0, u0]

r_start = 0.1
r_end = 100
r = np.linspace(r_start, r_end, 1000)

# Solve the differential equations
#sol = solve_ivp(dSdr, t_span=(min(r), max(r)), y0=S0, t_eval=r, method='BDF')
#sol = solve_ivp(dSdr, t_span=(min(r), max(r)), y0=S0, t_eval=r, method='RK45')
sol = solve_ivp(dSdr, t_span=(r_start, r_end), y0=S0, t_eval=r, method='LSODA', atol=1e-8, rtol=1e 8)



# Check if integration was successful
if sol.success:
    print("Integration successful")
else:
    print("Integration failed:", sol.message)

# Debugging information
print("Number of function evaluations:", sol.nfev)
print("Number of accepted steps:", sol.t_events)
print("Number of steps that led to errors:", sol.t)

# Plot results
plt.figure(figsize=(10, 6))
plt.plot(sol.t, sol.y[0], label='g(r)')
plt.plot(sol.t, sol.y[2], label='f(r)')
plt.xlabel('r')
plt.ylabel('Function values')
plt.legend()
plt.title('Solutions of the differential equations')
plt.show()

I would expect from the system that f goes exponentially towards 1 and g asymthotically towards 0 (with a little bumb in between), but this does not happen. Are there other methods that I can try? Is my code just wrong or am I missing something? It would be nice if someone could help me :)

Answer 1

You want to modify your initial conditions to avoid singularities. Make sure your equations are correct.

import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt


def dSdr(r, S):
    g, v, f, u = S
    dgdr = v
    dvdr = -2 / r * v + 2 / r ** 2 * g - 3 / r * e * g ** 2 - \
        e ** 2 * g ** 3 - e / (2 * r) * f**2 - e**2 / 4 * f ** 2 * g
    dfdr = u
    dudr = -2 / r * u + 2 / r ** 2 * f - 2 * e / r * f * g - e ** 2 / 2 * f * g**2 + mu ** 2 * f - lambd * f ** 3

    return [dgdr, dvdr, dfdr, dudr]


e = 1
mu = 1
lambd = 1

g0 = 0.0001
v0 = 0.0001
f0 = 1
u0 = 0.0001

S0 = [g0, v0, f0, u0]

r_start, r_end = 0.000001, 1
r = np.linspace(r_start, r_end, 1000)

sol = solve_ivp(dSdr, t_span=(r_start, r_end), y0=S0, t_eval=r, method='LSODA', atol=1e-8, rtol=1e-8)


# print(sol.t)


plt.figure(figsize=(10, 6))
plt.plot(sol.t, sol.y[0], label='g(r)')
plt.plot(sol.t, sol.y[2], label='f(r)')
plt.xlabel('r')
plt.ylabel('Function values')
plt.legend()
plt.title('Solutions of the differential equations')
plt.show()

Prints

---

To properly come up with a correct numerical solution, we do need to set all our boundary conditions and initial values correctly.

Here, instead of setting r = 0 and r = math.inf, we use small and large values for simulation:

r_start, r_end = 0.000001, 10000

Prints