pythonnumpynumericdifferential-equations

How to evaluating solutions to Coupled Nonlinear Elliptic ODEs with Relaxation Method?

I would like to determine the following system of coupled, nonlinear, elliptic ODEs of second order with boundary conditions f(0) = h(0) = 0, f(1) = h(1) = 1.

enter image description here

I used a relaxation method to solve the system (at least that's what I understand by it). However, my code gives me a rather caotic result, which cannot be correct. I am not used to dealing with such difficult differential equations, is another method perhaps more suitable? Perhaps someone sees a mistake that I have made.

import numpy as np
import matplotlib.pyplot as plt

# Define parameters
N = 100  # Number of points
lambda_val = 1  # Example value for lambda
g = 1  # Example value for g
tolerance = 1e-6  # Convergence tolerance
max_iterations = 10000  # Maximum number of iterations

# Discretize the domain
x = np.linspace(0, 1, N + 1)
dx = x[1] - x[0]

# Initial guesses
f = x.copy()
h = x.copy()

# Boundary conditions
f[0], h[0] = 0, 0
f[1], h[1] = 1, 1

# Relaxation method
for iteration in range(max_iterations):
    f_old = f.copy()
    h_old = h.copy()

    # Update f
    for i in range(1, N):
        # Discretized derivatives
        f_prime = (f[i + 1] - f[i - 1]) / (2 * dx)
        f_double_prime = (f[i + 1] - 2 * f[i] + f[i - 1]) / (dx ** 2)

        # Avoid division by zero
        h_term = x[i] ** 2 * h[i] ** 2 / (x[i] - 1) ** 2 if x[i] != 1 else 0
        rhs_f = (1 - f[i]) * (8 * (1 - 2 * f[i]) * f[i] - h_term) / 4
        lhs_f = 2 * (x[i] - 1) * x[i] ** 2 * f_prime + (x[i] - 1) ** 2 * x[i] ** 2 * f_double_prime

        # Ensure the update step is reasonable
        update_f = (rhs_f - lhs_f) * dx ** 2
        if np.abs(update_f) < 1e6:  # Avoid too large updates
            f[i] = f_old[i] + update_f

    # Update h
    for i in range(1, N):
        # Discretized derivatives
        h_prime = (h[i + 1] - h[i - 1]) / (2 * dx)
        h_double_prime = (h[i + 1] - 2 * h[i] + h[i - 1]) / (dx ** 2)

        # Avoid division by zero
        f_term = 2 * (f[i] - 1) ** 2
        lambda_term = x[i] ** 2 * lambda_val * (h[i] ** 2 - 1) / (g ** 2 * (x[i] - 1) ** 2) if x[i] != 1 else 0
        rhs_h = h[i] * (f_term + lambda_term)
        lhs_h = 2 * (x[i] - 1) ** 2 * x[i] * h_prime + (x[i] - 1) ** 2 * x[i] ** 2 * h_double_prime

        # Ensure the update step is reasonable
        update_h = (rhs_h - lhs_h) * dx ** 2
        if np.abs(update_h) < 1e6:  # Avoid too large updates
            h[i] = h_old[i] + update_h

    # Check for convergence
    if np.max(np.abs(f - f_old)) < tolerance and np.max(np.abs(h - h_old)) < tolerance:
        print(f"Converged after {iteration} iterations.")
        break
else:
    print("Did not converge within the maximum number of iterations.")

# Print the final solution
print("f =", f)
print("h =", h)

# Plotting the results
plt.figure(figsize=(12, 6))

# Plot f(x)
plt.subplot(1, 2, 1)
plt.plot(x, f, label='f(x)')
plt.xlabel('x')
plt.ylabel('f(x)')
plt.title('Solution for f(x)')
plt.legend()
plt.grid(True)

# Plot h(x)
plt.subplot(1, 2, 2)
plt.plot(x, h, label='h(x)', color='orange')
plt.xlabel('x')
plt.ylabel('h(x)')
plt.title('Solution for h(x)')
plt.legend()
plt.grid(True)

plt.tight_layout()
plt.show()

enter image description here

Solution

  • As I put in a comment, I don’t think you should do a transformation for x just to bring infinity to a finite point: I think you should stick with the original equations in r (equations 11.44 and 11.45 in your linked document) and just put the right boundary at a reasonably large number. I have solved these original equations by the finite-volume method in a slightly improved version of that for one of your previous threads, Solving coupled 2nd order ODE, numerical in Python

    Slightly rearranged – for reasons given below – the original equations are

    enter image description here

    These CONSERVATIVE equations can be discretised by integrating over a cell. Note that I divide the source term on the RHS into an explicit and an implicit part; the coefficient for the latter must be negative.

    enter image description here

    Each of these can then be written as a TRI-DIAGONAL SYSTEM:

    enter image description here

    where

    enter image description here

    Making the implicit source term negative ensures that this is diagonally dominant. The boundary conditions are usually handled through the source term.

    In this instance I found it best to solve the system iteratively using the TRI-DIAGONAL MATRIX ALGORITHM, a version of which is included in the code below.

    Code:

    import numpy as np
import matplotlib.pyplot as plt


#-----------------------------------------------------------


def tdma( a, b, c, d ):
#******************************************************************************
# Solve, using the TRI-DIAGONAL MATRIX ALGORITHM (TDMA)
#     a_i x_i-1 + b_i x_i + c_i x_i+1 = d_i,     i = 0, n-1
#
# Effectively, this is the n x n matrix equation, where
# a[i], b[i], c[i] are the non-zero diagonals of the matrix and d[i] is the rhs.
# a[0] and c[n-1] are not used.
#******************************************************************************
    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


#-----------------------------------------------------------


N = 400                                # number of cells (1-N); added boundaries 0 and N+1
rmin, rmax = 0.0, 20.0                 # minimum and maximum r
dr = rmax / N                          # cell width
fL, fR, hL, hR = 0, 1, 0, 1            # boundary conditions on f and h at L (left) and R (right) boundaries
lamda = 0.5


# Cell / node arrangement:
#
#       0                                          N+1
#       |  1     2                        N-1    N  |
#       |-----|-----| .  .  .  .  .   . |-----|-----|
#

r  = np.linspace( rmin - 0.5 * dr, rmax + 0.5 * dr, N + 2 )      # cell nodes (except for boundaries)
r[0] = rmin;   r[N+1] = rmax                                     # boundary nodes on faces of cells
rv = np.linspace( rmin, rmax, N + 1 )                            # cell faces (index is the right side of the cell)


# Connection coefficients
awf = np.zeros( N + 2 )
aef = np.zeros( N + 2 )
awh = np.zeros( N + 2 )
aeh = np.zeros( N + 2 )
# Main cells
aef[1:N+1] = 1              / dr ** 2
aeh[1:N+1] = rv[1:N+1] ** 2 / dr ** 2
awf[1:N+1] = aef[0:N]
awh[1:N+1] = aeh[0:N]
# Boundaries
awf[1], awh[1], aef[N], aeh[N] = 0, 0, 0, 0     # zero out the matrix coefficients
awfL = 2             / dr ** 2                  # apply boundary conditions via the source term
aefR = 2             / dr ** 2
awhL = 2 * rmin ** 2 / dr ** 2
aehR = 2 * rmax ** 2 / dr ** 2


# Initial values
f = np.zeros( N + 2 )
h = np.zeros( N + 2 )
# Boundary conditions
f[0] = fL;   f[N+1] = fR;   h[0] = hL;   h[N+1] = hR


alpha = 0.5                  # Under-relaxation factor
niter = 0
s  = np.zeros( N + 2 )
sp = np.zeros( N + 2 )
for _ in range( 2000 ):
    niter += 1

    ftarget = f.copy()
    htarget = h.copy()

    ######################
    # Update f equation  #
    ######################
    # Source term (write as S + Sp.f, where Sp is negative)
    for i in range( 1, N + 1):                             # cells 1 - N
       s [i] = + h[i] ** 2 / 4
       sp[i] = - h[i] ** 2 / 4
       if ( 1 - f[i] ) * ( 1 - 2 * f[i] ) > 0:
          sp[i] += -2        * ( 1 - f[i] ) * ( 1 - 2 * f[i] ) / r[i] ** 2
       else:
          s [i] += -2 * f[i] * ( 1 - f[i] ) * ( 1 - 2 * f[i] ) / r[i] ** 2
    # Boundary conditions (applied through the source term)
    s[1] += awfL * fL;   sp[1] -= awfL
    s[N] += aefR * fR;   sp[N] -= aefR
    # Diagonal term (including implicit part of the source)
    ap = awf + aef - sp
    # Solve tri-diagonal system
    ftarget[1:N+1] = tdma( -awf[1:N+1], ap[1:N+1], -aef[1:N+1], s[1:N+1] )


    ######################
    # Update h equation  #
    ######################
    # Source term (write as S + Sp.h, where Sp is negative)
    for i in range( 1, N + 1):                             # cells 1 - N
       s [i] =                             lamda * r[i] ** 2 * h[i]
       sp[i] = -( 2 * ( 1 - f[i] ) ** 2 +  lamda * r[i] ** 2 * h[i] ** 2 )
    # Boundary conditions (applied through the source term)
    s[1] += awhL * hL;   sp[1] -= awhL
    s[N] += aehR * hR;   sp[N] -= aehR
    ap = awh + aeh - sp
    # Solve tri-diagonal system
    htarget[1:N+1] = tdma( -awh[1:N+1], ap[1:N+1], -aeh[1:N+1], s[1:N+1] )


    change = ( np.linalg.norm( ftarget - f ) + np.linalg.norm( htarget - h ) ) / N

    # Under-relax the update
    f = alpha * ftarget + ( 1 - alpha ) * f
    h = alpha * htarget + ( 1 - alpha ) * h

    # print( niter, change )
    if change < 1.0e-10: break

print( niter, " iterations " )


plt.plot( r, f, label='f' )
plt.plot( r, h, label='h' )
plt.xlim( ( 0, 10  ) )
plt.ylim( ( 0, 1.1 ) )
plt.grid()
plt.legend()
plt.show()

    Output:

    enter image description here

    Compare Figure 11.6 in your linked document:

    enter image description here