BVP problem with conjugate gradient method

Question

Life and death situation. I have the differential equaction $$1.1 \dfrac{d^2T}{dz} - 0.01\dfrac{dT}{dz} = 0, T(-1000) = 0.12, T(0) = -1.55$$ that im trying to solve with conjugate gradient method.

So i translated the eq to the equation system $Ax = b$ approximating derivatives with central difference. This is my code and approach to solve it. Im pretty sure that the function for conjugate gradient is correct implemented (since i used it suceffully for another problem). I am aware that in order to use the method the matrix has to be positive definit which it is (i checked the eigen values). Anyway this is the code:

import numpy as np
import matplotlib.pyplot as plt

def createA(N):
    
    h= (1000)/ N  
    A = np.zeros((N+1,N+1)) #Kvadratisk n+1 matris med nollor
    A[0,0] = 1 #elementet på rad 1, kolumn 1
    A[N,N] = 1 # elementet på rad n+1, kolumn n+1
    p_1 = 1.1 + 0.005*h #Uppdaterar värdet på h^2 - 2 
    p_2 = 1.1 - 0.005*h
    
    for i in range(1, N):  #Sätter ut ettor på till vänster/höger om diagnoalelement på rad 2 till N
        A[i, i-1] = p_1
        A[i, i] = -2.2
        A[i, i+1] = p_2
        
    return(A)


def createF(N ,leftbc, rightbc):
    
    v = np.zeros(N+1) 
    v[0] = leftbc  
    v[N] = rightbc
    
    return(v)


def cg(A, b, x_cg, maxit):

    rk = b-np.dot(A,x_cg)
    p = rk
    
    for k in np.arange( maxit ):
        alpha = np.dot(rk,rk)/np.dot(p , np.dot(A , p ) )
        x_cg = x_cg + alpha *p
        rkp1 = rk - alpha*np.dot(A,p)
        beta = np.dot(rkp1, rkp1)/np.dot(rk,rk)
        p = rkp1 + beta*p
        rk = rkp1
    
    return x_cg

N = 5
matrix = createA(N)
vektor = createF(N,0.12,-1.55)
guess = np.zeros(N+1)

solution = cg(matrix, vektor, guess,100)

but this just dont work, i have googled the whole internet and cant find the solution. When i change the number of iterations or N the solution just goes crazy. Maybe there is a smart guy that can help.

Just in case: This is the analytical solution for the eq: $$ T(z) = -1.55-\dfrac{1.67}{\exp\left(\dfrac{-10}{1.1}\right)-1} + \dfrac{1.67}{\exp\left(\dfrac{-10}{1.1}\right)-1} \exp\left(\dfrac{0.01z}{1.1}\right) $$

Thanks in advance fellows! #Time_is_a_factor

Answer 1

I'm almost sure you can use the conjugate gradient method but it would be easier to convert the differential equation into a boundary value problem.

Something like this

import numpy as np
from scipy.integrate import solve_bvp
import matplotlib.pyplot as plt

# Define the differential equation and boundary conditions
def fun(z, T):
    dTdz = T[1]
    d2Tdz2 = 0.01/1.1 * dTdz 
    return np.vstack((dTdz, d2Tdz2))

def bc(Ta, Tb):
    return np.array([Ta[0] - 0.12, Tb[0] + 1.55])

# Define the mesh for the variable (z)
z_mesh = np.linspace(-1000, 0, 1000)

# Initial guess for the solution
initial_guess = np.zeros((2, z_mesh.size))

# Solve the boundary value problem using solve_bvp
solution = solve_bvp(fun, bc, z_mesh, initial_guess)

# Extract the solution
T_solution = solution.y[0]

def analytical(z):
    return -1.55-1.67/(np.exp(-10/1.1)-1)+1.67/(np.exp(-10/(1.1))-1)*np.exp(0.01*z/(1.1))

print(T_solution[-1])
# Plot the solution
plt.figure(figsize=(8, 6))
plt.plot(z_mesh,analytical(z_mesh),label="Analytical solution")

plt.plot(z_mesh, T_solution, '--',label='Numerical solution')
plt.xlabel('z')
plt.ylabel('T(z)')
plt.title('Solution of the Differential Equation')
plt.legend()
plt.grid(True)
plt.show()

As you can see this matches the analytical solution.