How to solve three coupled differential equations in python using RK-4 and shooting method? or using solve_bvp?
I am trying to solve three coupled differential equations in Python. I am using RK-4 techniques with Shooting method. I am trying to plot the f and N functions.
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import fsolve
alpha = 0.1
def F(r, f, v, N, sig):
g = np.sin(f)
h = np.cos(f)
i = np.sin(2 * f)
dfdr = v
dNdr = -(alpha*r**2 * N * v**2 *(r**2 + 2 * g**2)+alpha * g**2 * (g**2 +2* r**2) + r**2*(N-1))/r**3
dsigdr = ( alpha * sig * v**2 * (r**2 + 2 * g**2)) / r
f1= v*((2 * dsigdr * g**2)/(sig * r**2)+(2 * dNdr * g**2)/(N * r**2)+((dsigdr)/sig +dNdr / N + 2 / r))
dvdr = -(f1 -i/(N * r**2)+(-(2 * g**3 * h)/(N * r**4)+(v**2 * i)/r**2)) / (1 + (2 * g**2) / r**2)
return [dfdr, dNdr, dsigdr, dvdr]
def rk4(F, r_range, f0, v0, N0, s0, h):
r_values = np.arange(r_range[0], r_range[1], h)
f_values = np.linspace(f0, 0, len(r_values))
v_values = np.linspace(v0, 0, len(r_values))
N_values = np.linspace(N0, 1, len(r_values))
S_values = np.linspace(s0, 1, len(r_values))
for i in range(1, len(r_values)):
r = r_values[i-1]
f = f_values[i-1]
v = v_values[i-1]
N = N_values[i-1]
sig = S_values[i-1]
k1f, k1v, k1N, k1S = F(r, f, v, N, sig)
k2f, k2v, k2N, k2S = F(r + 0.5*h, f + 0.5*k1f*h, v + 0.5*k1v*h, N + 0.5*k1N*h, sig +0.5*k1S*h)
k3f, k3v, k3N, k3S = F(r + 0.5*h, f + 0.5*k2f*h, v + 0.5*k2v*h, N + 0.5*k2N*h, sig +0.5*k2S*h)
k4f, k4v, k4N, k4S = F(r + h, f + k3f*h, v + k3v*h, N + k3N*h, sig + k3S*h)
f_values[i] = f + ((k1f + 2*k2f + 2*k3f + k4f)*h) / 6
v_values[i] = v + ((k1v + 2*k2v + 2*k3v + k4v)*h) / 6
N_values[i] = N + ((k1N + 2*k2N + 2*k3N + k4N)*h) / 6
S_values[i] = sig + ((k1S + 2*k2S + 2*k3S + k4S)*h)/ 6
return r_values, f_values, v_values, N_values, S_values
rmin = 0.001
rmax = 4
r_range=[rmin, rmax]
h = 0.00001
def objective(u):
v0, s0 = u
r_values, f_values, v_values, N_values, S_values = rk4(F, r_range, np.pi, v0, 1, s0, h)
y1 = f_values[-1]
y2 = S_values[-1]
print(f"Objective results: y1={y1}, y2={y2}")
return [y1 , y2 - 1]
initial_guess=[-0.478, -1.216]
v1 = fsolve(objective, initial_guess, xtol=1e-8)
print("The initial values of v0 and sig0 given by fsolve:", v1[0], " and ", v1[1])
r_values, f_values, v_values, N_values, S_values = rk4(F, [rmin, rmax], np.pi, v1[0], 1, v1[1], h)
y3 = f_values[-1]
y4 = S_values[-1]
print("The value of f and sig at outer boundary is:", y3, y4)
if np.abs(y3) < 1e-6 and 1<= y4 <1.1:
print("Boundary condition f(rmax) = 0 is satisfied, and the value of f at outer boundar is:", y3)
print("Boundary condition sig(rmax) = 0 is satisfied, and the value of sig at outer boundar is:", y4)
else:
print("The solutions did not converge")
plt.figure(figsize=(14, 10))
plt.plot(r_values, f_values, label='f(r)')
plt.plot(r_values, N_values, label='N(r)')
plt.title(r"$\alpha=0.1$")
plt.xlabel('r')
plt.grid()
plt.legend()
plt.show()
But the nature of the plot is completely different from my expectations. The output is
Objective results: y1=0.000372840219254553, y2=1.0001907509516712
Objective results: y1=0.000372840219254553, y2=1.0001907509516712
Objective results: y1=0.000372840219254553, y2=1.0001907509516712
Objective results: y1=0.00037294172854991653, y2=1.0001906882951273
Objective results: y1=0.000372824083910513, y2=1.0001907727358128
Objective results: y1=1.8398001535383646e-07, y2=0.9999999172810013
Objective results: y1=1.1072732970996606e-11, y2=1.0000000000024691
Objective results: y1=2.1284157569008075e-12, y2=0.9999999999969905
Objective results: y1=1.1784771089214975e-13, y2=0.9999999999997673
The initial values of v0 and sig0 given by fsolve: -0.4780664768861783 and -1.2166450055671914
The value of f and sig at the outer boundary is: 1.1784771089214975e-13 0.9999999999997673
The solutions did not converge
As you see the value of sig at the outer boundary is almost one (0.9999999999997673).
These are the three coupled differential equations with boundary conditions
I forgot to mention an extra condition on function N, here an extra condition on function N is given
I also tried the above code with a smaller value of h (=0.000001) and different values of rmin, rmax, alpha and initial_guess but the result was the same.
I also tried to solve the problem using solve_bvp. Here is the code using solve_bvp
import numpy as np
from scipy.integrate import solve_bvp
import matplotlib.pyplot as plt
alpha = 0.1
def F(r, f):
#f, v, N, sig = S
g = np.sin(f[0])
h = np.cos(f[0])
i = np.sin(2 * f[0])
#dfdr = v
dNdr = -( alpha * r**2 * f[2] * f[1]**2 * (r**2 + 2* g**2) + alpha * g**2 * (g**2 +2* r**2) + r**2 * (f[2] - 1)) / r**3
dsigdr = ( alpha * f[3] * f[1]**2 * (r**2 + 2 * g**2)) / r
ter1=((dsigdr) / f[3] + dNdr / f[2] + 2 / r +(2 * dsigdr * g**2) / (f[3] * r**2) + (2 * dNdr* g**2) / (f[2] * r**2))
dvdr = -((f[1] * ter1 - i / (f[2] * r**2)) + (-(2 * g**3 * h) / (f[2] * r**4) + (f[1]**2 * i) / r**2)) / (1 + (2 * g**2) / r**2)
return np.vstack((f[1], dvdr, dNdr, dsigdr))
def bc(fa, fb):
return np.array([fa[0]-np.pi, fb[0]-0, fa[2]-1, fb[3]-1])
rmin = 0.001
rmax = 5
x = np.linspace(rmin, rmax, 100000)
y=np.ones([4, x.size])
y[0, :] = np.linspace(np.pi, 0, x.size)
y[1, :] = np.linspace(-0.2, 0, x.size)
y[2, :] = np.linspace(1, 1, x.size)
y[3] = np.ones(x.size)
sol = solve_bvp(F, bc, x, y, tol=1e-8, max_nodes=100000, verbose=2)
if sol.success:
print("Solution converged")
else:
print("Solution did not converge")
print("Residuals:", sol.message)
print("Residuals max:", np.max(sol.rms_residuals))
print(sol)
#Plot the solutions
plt.figure(figsize=(14, 10))
plt.plot(sol.x, sol.y[0], label='f(r)')
#plt.plot(sol.x, sol.y[1], label='V(r)')
plt.plot(sol.x, sol.y[2], label='N(r)')
plt.xlabel('r')
plt.grid()
plt.legend()
plt.show()
When I am using solve_bvp, I am getting the following error
Iteration Max residual Max BC residual Total nodes Nodes added
1 1.06e+00 3.41e-21 100000 (199998)
Number of nodes is exceeded after iteration 1.
Maximum relative residual: 1.06e+00
Maximum boundary residual: 3.41e-21
Solution did not converge
Residuals: The maximum number of mesh nodes is exceeded.
Residuals max: 1.0570652692628755
message: 'The maximum number of mesh nodes is exceeded.'
niter: 1
p: None
rms_residuals: array([4.25858920e-01, 2.54417888e-01, 1.98041640e-01, ...,
6.12622641e-05, 6.11022248e-05, 6.09430998e-05])
sol: <scipy.interpolate._interpolate.PPoly object at 0x000001B21E651DF0>
status: 1
success: False
The plot using solve_bvp is better but the success status is false.
Solution
I am referring to your "shooting" code.
You have a major numerical error in the derivative routine, where return [dfdr, dNdr, dsigdr, dvdr] is the wrong order. The return order should be
return [dfdr, dvdr, dNdr, dsigdr]
You have an error in your equation (1), where you are missing a factor of 2 multiplying sin(2f) in the second line. Referring to the second of the references that you quote, this comes from differentiating 2.sin^2(f), which gives 4.sin(f)cos(f).df/dr, or 2.sin(2f).df/dr.
Your whole derivative routine is difficult to align with the original equations, so I have taken the liberty of rewriting it below.
You have correctly moved the lower boundary away from r=0. However, you can use the Taylor-series/power-series expansion to improve the values at rmin in this case: f(rmin)=pi+(df/dr)*rmin.
I think that you should increase rmax (to 10, say), whilst with Runge-Kutta you can easily increase h by at least a factor of 10 to speed up results.
Complete code below. Note that I have reverted to a value alpha=0.04 so that you can compare directly with your expected graph. (I have also divided f by pi in the picture, as is done there.)
Final code:
import math
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import fsolve
alpha = 0.04
def F(r, f, v, N, sig):
S = math.sin(f)
S2 = math.sin(2 * f)
dfdr = v
dNdr = -alpha*N*(r**2+2*S**2)*v**2/r - alpha*S**2*(S**2+2*r**2)/r**3 - (N-1)/r
dsigdr = alpha*sig*v**2*(r**2+2*S**2)/r
dvdr = ( (dNdr/N + dsigdr/sig)*(1+2*S**2/r**2) + 2/r + 2*S2*v/r**2 ) * v -(S2/N/r**2)*(1+N*v**2+S**2/r**2)
dvdr = -dvdr/(1+2*S**2/r**2)
return [dfdr, dvdr, dNdr, dsigdr]
def rk4(F, r_range, f0, v0, N0, s0, h):
r_values = np.arange(r_range[0], r_range[1], h)
f_values = np.linspace(f0, 0, len(r_values))
v_values = np.linspace(v0, 0, len(r_values))
N_values = np.linspace(N0, 1, len(r_values))
S_values = np.linspace(s0, 1, len(r_values))
for i in range(1, len(r_values)):
r = r_values[i-1]
f = f_values[i-1]
v = v_values[i-1]
N = N_values[i-1]
sig = S_values[i-1]
k1f, k1v, k1N, k1S = F(r, f, v, N, sig)
k2f, k2v, k2N, k2S = F(r + 0.5*h, f + 0.5*k1f*h, v + 0.5*k1v*h, N + 0.5*k1N*h, sig +0.5*k1S*h)
k3f, k3v, k3N, k3S = F(r + 0.5*h, f + 0.5*k2f*h, v + 0.5*k2v*h, N + 0.5*k2N*h, sig +0.5*k2S*h)
k4f, k4v, k4N, k4S = F(r + h, f + k3f*h, v + k3v*h, N + k3N*h, sig + k3S*h)
f_values[i] = f + ((k1f + 2*k2f + 2*k3f + k4f)*h) / 6
v_values[i] = v + ((k1v + 2*k2v + 2*k3v + k4v)*h) / 6
N_values[i] = N + ((k1N + 2*k2N + 2*k3N + k4N)*h) / 6
S_values[i] = sig + ((k1S + 2*k2S + 2*k3S + k4S)*h)/ 6
return r_values, f_values, v_values, N_values, S_values
rmin = 0.0001
rmax = 10
r_range=[rmin, rmax]
h = 0.0001
def objective(u):
v0, s0 = u
r_values, f_values, v_values, N_values, S_values = rk4(F, r_range, np.pi + v0 * rmin, v0, 1, s0, h)
y1 = f_values[-1]
y2 = S_values[-1]
print(f"Objective results: y1={y1}, y2={y2}")
return [y1 , y2 - 1]
initial_guess=[-2.794, 0.471 ]
v1 = fsolve(objective, initial_guess, xtol=1e-8)
print("The initial values of v0 and sig0 given by fsolve:", v1[0], " and ", v1[1])
r_values, f_values, v_values, N_values, S_values = rk4(F, [rmin, rmax], np.pi + v1[0] * rmin, v1[0], 1, v1[1], h)
y3 = f_values[-1]
y4 = S_values[-1]
print("The value of f and sig at outer boundary is:", y3, y4)
if abs(y3) + abs( y4 - 1 ) < 1e-6:
print("Boundary condition f(rmax) = 0 is satisfied, and the value of f at outer boundary is:", y3)
print("Boundary condition sig(rmax) = 0 is satisfied, and the value of sig at outer boundary is:", y4)
else:
print("The solutions did not converge")
plt.figure(figsize=(8, 6))
plt.plot(r_values, f_values/np.pi, label='f(r)/pi')
plt.plot(r_values, N_values, label='N(r)')
plt.title(r"$\alpha=0.04$")
plt.xlabel('r')
plt.grid()
plt.legend()
plt.show()
Output:
Objective results: y1=0.00019542960351076646, y2=0.9993130868376403
Objective results: y1=0.00019542960351076646, y2=0.9993130868376403
Objective results: y1=0.00019542960351076646, y2=0.9993130868376403
Objective results: y1=0.0001954504980834034, y2=0.9993130650610034
Objective results: y1=0.00019542960351076646, y2=0.9993131017285548
Objective results: y1=2.1493424939744863e-07, y2=1.000000127532956
Objective results: y1=5.273386299353826e-10, y2=1.0000000001353808
Objective results: y1=7.569440649813437e-13, y2=0.9999999999992873
The initial values of v0 and sig0 given by fsolve: -2.7943898358693358 and 0.47122759362590755
The value of f and sig at outer boundary is: 7.569440649813437e-13 0.9999999999992873
Boundary condition f(rmax) = 0 is satisfied, and the value of f at outer boundary is: 7.569440649813437e-13
Boundary condition sig(rmax) = 0 is satisfied, and the value of sig at outer boundary is: 0.9999999999992873
- AttributeError: install_layout when attempting to install a package in a virtual environment
- Python list comprehension - want to avoid repeated evaluation
- Hash algorithm for dynamic growing/streaming data?
- matplotlib - making labels for violin plots
- Python How to I check if last element has been reached in iterator tool chain?
- Polars and the Lazy API: How to drop columns that contain only null values?
- Why are my Mean, Var, and Std outputs from NumPy different from what the online grader expects?
- Correlation dataframe convertion from results from pl.corr
- Polars DataFrame transformation
- Discord rate limiting while only sending 1 request per minute
- Check if column contains (/,-,_, *or~) and split in another column - Pandas
- How to draw a rectangle at (x,y) in a PyQt GraphicsView?
- how to calculate correlation between ten columns with polars
- How to set class attribute with await in __init__
- Detect hindi encoding, response received from Facebook API in Python
- Is it possible to write a horizontal if statement with a multi-line body?
- Max length of items in list
- Cannot subclass multiprocessing Queue in Python 3.5
- How can I get notified of updates to Python packages in a unified way?
- Using python AST to traverse code and extract return statements
- merge groups of columns in a polars dataframe to single columns
- Group Pandas DataFrame by Continuous Date Ranges
- Flask login @login_required not working
- Odoo: one2many and many2one? KeyError:'___'
- merge some columns in a Polars dataframe and duplicate the others
- Python: Create table from string mixed with separators using FOR loops
- How do I type hint a method with the type of the enclosing class?
- How can I verify an emails DKIM signature in Python?
- Writing a class that accepts a callback in Python?
- Python Paramiko channel.exec_command not returning output intermittently