Solving non-linear ODE system with boundary conditions in Python

Question

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

If H=grF as you have written then the last bracket of the H equation would be identically 0. Is this what you intended? — lastchance
– lastchance, Commented Oct 16, 2024 at 12:47
@lastchance Ohh sorry for my bad notation. F(r) is a function and F is a constant. F does not appear in the code, I have set F = 1, only the function F(r) appears. The term is therefore not zero. — Hendriksdf5
– Hendriksdf5, Commented Oct 16, 2024 at 16:36
Are you (@Hendriksdf5) trying to solve equations for K and H or for W and F? Your code doesn't match the equations listed. Also, what are the intended boundary conditions at r = 0? — lastchance
– lastchance, Commented Oct 17, 2024 at 19:52

lastchance · Accepted Answer · 2024-10-17 20:30:18Z

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:

Thank you very much! The method of calculating the system in the original functions and then converting them seems to make more sense to me! Do you happen to know a good textbook that explains how to solve such coupled second order systems? I come across such systems relatively often and it would be nice if I could do it myself.
I am coding it from the way that I would approach it via the finite-volume method for CFD (although there is only a diffusion/viscous term here, not advection). Thus a simple book on CFD via the finite-volume method might help you (indirectly).

Collectives™ on Stack Overflow

Solving non-linear ODE system with boundary conditions in Python

1 Answer 1

2 Comments

Your Answer

Hot Network Questions

Collectives™ on Stack Overflow

1 Answer 1

2 Comments

Your Answer

Sign up or log in

Post as a guest

Related