# BVP problem with conjugate gradient method

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

• Unless you have some modifications, CG requires that $A$ be positive definite (symmetric). Commented Sep 19, 2023 at 1:59
• Yes, i wrote that te matrix is positive definit Commented Sep 19, 2023 at 8:41

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.