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

• Hendriksdf5: dSdr_scaled describes the same equation only scaled to some factor, I thought if I get the constants out it could be simpler
• Robert Dodier: Well, are the derivatives always zero for u = v = f = g = 0 ? If so, you could keep the initial conditions and consider if maybe there should be an additional inhomogeneous term, or keep the equations and try nonzero initial conditions. You might consider trying math.stackexchange.com instead since it seems like this problem is more about concepts than programming.
• lastchance: Can you show a link to where your equations come from - there's quite a lot of discrepancy between your post and your code.

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

enter image description here

---

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

enter image description here

• Hendriksdf5: Thanks for the effort, at least there is now a solution, but it doesn't look like it should. My boundary con. should be f(0) = g(0) = 0 and f(infinity) = mu/sqrt(lambd) g(infinity) = 0 so that the system makes physical sense. Do you happen to know how I can imply something like this in the code or what I have to look for to learn it?
• Aicody: @Hendriksdf5 Updated. Make sure your equations, numbers and boundaries are all correct.
• lastchance: @Hendriksdf5, the boundary conditions you have just stated make it a boundary-value problem, not an initial-value problem. Please provide a link to the original equations.