Solving non-linear ODE system with boundary conditions in Python
I am trying to solve this system of differential equations:
The functions should behave like this, as r-> infinity (F = g = lambda = const.):
A while ago I had a similar problem and a user was kind enough to provide me with an algorithm with which I could solve this ODS. Now I have tried to apply the TDM algorithm to the new problem but I only get the trivial solution. Does anyone know what I have done wrong or how I can solve this problem?
import numpy as np
import matplotlib.pyplot as plt
# TDMA function to solve the tri-diagonal matrix equation
def tdma(a, b, c, d):
n = len(d)
p = np.zeros(n)
q = np.zeros(n)
x = np.zeros(n)
# Forward pass
i = 0
denominator = b[i]
p[i] = -c[i] / denominator
q[i] = d[i] / denominator
for i in range(1, n):
denominator = b[i] + a[i] * p[i - 1]
if abs(denominator) < 1.0e-10:
print("No solution")
return x
p[i] = -c[i] / denominator
q[i] = (d[i] - a[i] * q[i - 1]) / denominator
# Backward pass
x[n - 1] = q[n - 1]
for i in range(n - 2, -1, -1):
x[i] = p[i] * x[i + 1] + q[i]
return x
# Parameters
N = 1000 # Number of cells (1-N); added boundaries 0 and N+1
rmin, rmax = 0.0, 100.0 # Minimum and maximum r
dr = rmax / N # Cell width
lamda = 10
g = 1.0 # You can adjust this parameter if needed
# Boundary conditions on F and W
FL, FR = 0, 1
WL, WR = 0, 1 / rmax
# Cell / node arrangement
r = np.linspace(rmin - 0.5 * dr, rmax + 0.5 * dr, N + 2)
r[0], r[N + 1] = rmin, rmax
rv = np.linspace(rmin, rmax, N + 1)
# Coefficients
awF = np.zeros(N + 2)
aeF = np.zeros(N + 2)
awW = np.zeros(N + 2)
aeW = np.zeros(N + 2)
# Main cells
aeF[1:N + 1] = 1 / dr ** 2
aeW[1:N + 1] = rv[1:N + 1] ** 2 / dr ** 2
awF[1:N + 1] = aeF[0:N]
awW[1:N + 1] = aeW[0:N]
# Boundaries
awF[1], awW[1], aeF[N], aeW[N] = 0, 0, 0, 0
awFL = 2 / dr ** 2
aeFR = 2 / dr ** 2
awWL = 2 * rmin ** 2 / dr ** 2
aeWR = 2 * rmax ** 2 / dr ** 2
# Initial values
F = np.ones(N + 2) * 1e-5
W = np.ones(N + 2) * 1e-5
# Boundary conditions
F[0], F[N + 1] = FL, FR
W[0], W[N + 1] = WL, WR
# Under-relaxation factor
alpha = 0.5
niter = 0
s = np.zeros(N + 2)
sp = np.zeros(N + 2)
# Iteration loop
for _ in range(2000):
niter += 1
F_target = F.copy()
W_target = W.copy()
######################
# Update W equation #
######################
for i in range(1, N + 1):
K = 1 - g * r[i] * W[i]
H = g * r[i] * F[i]
# Clip values to prevent overflow
K = np.clip(K, -1e5, 1e5)
H = np.clip(H, -1e5, 1e5)
s[i] = K * (K**2 - 1) + H**2 * K
sp[i] = -2 * (1 - g * r[i]) / r[i]**2
if np.isnan(K) or np.isnan(H):
print(f"Numerical issue at r = {r[i]}, K = {K}, H = {H}")
break
s[1] += awFL * FL
sp[1] -= awFL
s[N] += aeFR * FR
sp[N] -= aeFR
ap = awF + aeF - sp
F_target[1:N + 1] = tdma(-awF[1:N + 1], ap[1:N + 1], -aeF[1:N + 1], s[1:N + 1])
######################
# Update F equation #
######################
for i in range(1, N + 1):
K = 1 - g * r[i] * W[i]
H = g * r[i] * F[i]
# Clip values to prevent overflow
K = np.clip(K, -1e5, 1e5)
H = np.clip(H, -1e5, 1e5)
s[i] = 2 * H * K**2 + lamda * H * ((H**2 / g**2) - r[i]**2 * F[i]**2)
if np.isnan(K) or np.isnan(H):
print(f"Numerical issue at r = {r[i]}, K = {K}, H = {H}")
break
s[1] += awWL * WL
sp[1] -= awWL
s[N] += aeWR * WR
sp[N] -= aeWR
ap = awW + aeW - sp
W_target[1:N + 1] = tdma(-awW[1:N + 1], ap[1:N + 1], -aeW[1:N + 1], s[1:N + 1])
# Calculate change
change = (np.linalg.norm(F_target - F) + np.linalg.norm(W_target - W)) / N
# Under-relax the update
F = alpha * F_target + (1 - alpha) * F
W = alpha * W_target + (1 - alpha) * W
if change < 1.0e-10:
break
print(niter, " iterations ")
# Plot the results
plt.plot(r, F, label='F(r)', color='g')
plt.plot(r, W, label='W(r)', color='r')
plt.legend(loc="upper right")
plt.xlabel("r")
plt.ylim(0)
plt.xlim(0, 20)
plt.grid(True)
plt.show()
Solution
I think that you should solve for W and F rather than H and K. If you transform the equations and rearrange you get (I think, but you had better check my algebra):
These are each discretised conservatively in the finite-volume form:
(where, for negative feedback, sp must be negative). This is then written in standard form as a tridiagonal system:
where
This can then be solved by the tridiagonal matrix algorithm (TDMA). The whole must be iterated because (a) the equations are coupled; and (b) the coefficients vary with the solution.
Code:
import numpy as np
import matplotlib.pyplot as plt
#-----------------------------------------------------------------------
# TDMA function to solve the tri-diagonal matrix equation
def tdma(a, b, c, d):
n = len(d)
p = np.zeros(n)
q = np.zeros(n)
x = np.zeros(n)
# Forward pass
i = 0
denominator = b[i]
p[i] = -c[i] / denominator
q[i] = d[i] / denominator
for i in range(1, n):
denominator = b[i] + a[i] * p[i - 1]
if abs(denominator) < 1.0e-10:
print("No solution")
return x
p[i] = -c[i] / denominator
q[i] = (d[i] - a[i] * q[i - 1]) / denominator
# Backward pass
x[n - 1] = q[n - 1]
for i in range(n - 2, -1, -1):
x[i] = p[i] * x[i + 1] + q[i]
return x
#-----------------------------------------------------------------------
# Parameters
N = 10000 # Number of cells (1-N); added boundaries 0 and N+1
rmin, rmax = 0.0, 100.0 # Minimum and maximum r
dr = rmax / N # Cell width
lamda = 10
g = 1.0 # You can adjust this parameter if needed
f = 1.0 # The "constant" value of F
# Boundary conditions on W and F
WL, WR = 0, 0
FL, FR = 0, f
# Cell / node arrangement
r = np.linspace(rmin - 0.5 * dr, rmax + 0.5 * dr, N + 2); r[0], r[N + 1] = rmin, rmax
rv = np.linspace(rmin, rmax, N + 1)
# Matrix coefficients (West and East, for each variable)
awW = np.zeros(N + 2)
aeW = np.zeros(N + 2)
awF = np.zeros(N + 2)
aeF = np.zeros(N + 2)
# Main cells
aeW[1:N + 1] = rv[1:N + 1] ** 2 / dr ** 2
aeF[1:N + 1] = rv[1:N + 1] ** 2 / dr ** 2
awW[1:N + 1] = aeW[0:N]
awF[1:N + 1] = aeF[0:N]
# Boundaries
awW[1], awF[1], aeW[N], aeF[N] = 0, 0, 0, 0
awWL = 0
aeWR = 2 * rmax ** 2 / dr ** 2
awFL = 0
aeFR = 2 * rmax ** 2 / dr ** 2
# Initial values
W = np.ones(N + 2) * 0.5
F = np.ones(N + 2) * 0.5
# Boundary conditions
W[0], W[N + 1] = WL, WR
F[0], F[N + 1] = FL, FR
# Under-relaxation factor
alpha = 0.5
niter = 0
s = np.zeros(N + 2)
sp = np.zeros(N + 2)
# Iteration loop
for _ in range(20000):
niter += 1
W_target = W.copy()
F_target = F.copy()
######################
# Update W equation #
######################
s = g * r * ( 3 * W ** 2 + F ** 2 )
sp = -( 2 + g ** 2 * r ** 2 * ( W ** 2 + F ** 2 ) )
s[1] += awWL * WL; sp[1] -= awWL
s[N] += aeWR * WR; sp[N] -= aeWR
ap = awW + aeW - sp
W_target[1:N + 1] = tdma(-awW[1:N + 1], ap[1:N + 1], -aeW[1:N + 1], s[1:N + 1])
######################
# Update F equation #
######################
s = lamda * r ** 2 * f ** 2 * F
sp = -( 2 * ( 1 - g * r * W ) ** 2 + lamda * r ** 2 * F ** 2 )
s[1] += awFL * FL; sp[1] -= awFL
s[N] += aeFR * FR; sp[N] -= aeFR
ap = awF + aeF - sp
F_target[1:N + 1] = tdma(-awF[1:N + 1], ap[1:N + 1], -aeF[1:N + 1], s[1:N + 1])
# Calculate change
change = (np.linalg.norm( W_target - W ) + np.linalg.norm( F_target - F ) ) / N
# Under-relax the update
W = alpha * W_target + ( 1 - alpha ) * W
F = alpha * F_target + ( 1 - alpha ) * F
if change < 1.0e-10: break
print(niter, " iterations ")
# Plot the results
plt.plot(r, F, label='F(r)', color='g')
plt.plot(r, W, label='W(r)', color='r')
plt.legend(loc="upper right")
plt.xlabel("r")
plt.ylim(0, 1.5)
plt.xlim(0, 10)
plt.grid(True)
plt.show()
Solution:
- how to average in a specific dimension with numpy.mean?
- Removed Realm, but still getting this error: module importing failed: invalid token (rlm_lldb.py, line 37) File "temp.py", line 1,
- Bool object does not support item assignment
- Conda env vs venv / pyenv / virtualenv / poetry / docker, etc
- How to import values from excel with pandas into tkinter faster?
- Speed up concatenation of excel files with Pandas
- Faster way to read Excel files to pandas dataframe
- Create Azure TimerTrigger Durable Function in python
- Read same file from src and test in PyCharm
- How to use Cython to compile Python 3 into C
- Outputting variable from one thread to other
- How to POST a JSON having a single body parameter in FastAPI?
- touchscreen raspberry pi phantom touch (reset fixes)
- How to determine the encoding of text
- Specify extras_require with pip install -e
- Plotting a second scaled y axis in matplotlib from one set of data
- How can I call scikit-learn classifiers from Java?
- How does Pandas.Series.nbytes work for strings? Results don't seem to match expectations
- How to install ursina on pydroid?
- auditlog with Django and DRF
- Python: How to RECURSIVELY remove None values from a NESTED data structure (lists and dictionaries)?
- How to Mark Repeated Entries as True Starting from the Second Occurrence Using NumPy?
- Python3.5 send object over socket (Pygame cam image)
- Value based partial slicing with non-existing keys is now deprecated
- Camelot: DeprecationError: PdfFileReader is deprecated
- How can I partially fill a numpy array with a value given a range?
- Conda wont install pdfplumber
- tomli-w, adding arrays from a python script
- Override value in pydantic model with environment variable
- Polars-Python: Compare list columns