Solving a tridiagonal matrix in python

Question

I have been looking at numerical methods to solve differential equations for chemical reactions. Usually I put the differential equation into a tridiagonal matrix using finite difference method, and then using a column vector for the boundary conditions. Once I have the matrix and vector I use scipy's linalg. However after formulating the tridiagonal matrix above I have no idea how to solve this on python since now the algebraic values are in the tridiagonal matrix, do I use some iterative method?

Any guidance will be greatly appreciated.

Answer 1

So I decided to use newton method for a system of equations to solve this problem, as recommended by @LutzLehmann.'J' is the Jacobian matrix and f is original matrix. This is not very efficient code but it got the job done.

guess = np.array([4,4,4,4,4,4,4,4,4,4,4])
for i in range(10):
     J = np.array([[-20003.0002-0.2*guess[0],1.0002,0,0,0,0,0,0,0,0,0],
              [1.0001,-1.0002-0.2*guess[1],0.0001,0,0,0,0,0,0,0,0],
              [0,1.0001,-1.0002-0.2*guess[2],0.0001,0,0,0,0,0,0,0],
              [0,0,1.0001,-1.0002-0.2*guess[3],0.0001,0,0,0,0,0,0],
              [0,0,0,1.0001,-1.0002-0.2*guess[4],0.0001,0,0,0,0,0],
              [0,0,0,0,1.0001,-1.0002-0.2*guess[5],0.0001,0,0,0,0],
              [0,0,0,0,0,1.0001,-1.0002-0.2*guess[6],0.0001,0,0,0],
              [0,0,0,0,0,0,1.0001,-1.0002-0.2*guess[7],0.0001,0,0],
              [0,0,0,0,0,0,0,1.0001,-1.0002-0.2*guess[8],0.0001,0],
              [0,0,0,0,0,0,0,0,1.0001,-1.0002-0.2*guess[9],0.0001],
              [0,0,0,0,0,0,0,0,0,1.0002,-1.0002-0.2*guess[10]]])
              f1=-20003.0002*guess[0]-0.1*(guess[0]**2) + 1.0002*guess[1]+20002
              f2= 1.0001*guess[0] -1.0002*guess[1]-0.1* (guess[1]**2)+0.0001*guess[2]
              f3 = 1.0001*guess[1]-1.0002*guess[2]-0.1*(guess[2]**2) +0.0001*guess[3]
              f4 = 1.0001*guess[2]-1.0002*guess[3]-0.1*(guess[3]**2) +0.0001*guess[4]
              f5 = 1.0001*guess[3]-1.0002*guess[4]-0.1*(guess[4]**2) +0.0001*guess[5]
              f6 = 1.0001*guess[4]-1.0002*guess[5]-0.1*(guess[5]**2) +0.0001*guess[6]
              f7 = 1.0001*guess[5]-1.0002*guess[6]-0.1*(guess[6]**2) +0.0001*guess[7]
              f8 = 1.0001*guess[6]-1.0002*guess[7]-0.1*(guess[7]**2) +0.0001*guess[8]
              f9 = 1.0001*guess[7]-1.0002*guess[8]-0.1*(guess[8]**2) +0.0001*guess[9]
             f10 =1.0001*guess[8]-1.0002*guess[9]-0.1*(guess[9]**2) +0.0001*guess[10]
             f11 = 1.0002*guess[9]-1.0002*guess[10]-0.1*(guess[10]**2) 
             f = np.array([f1,f2,f3,f4,f5,f6,f7,f8,f9,f10,f11])
             delta = linalg.solve(J,-f)
             guess = delta + guess
guess
array([0.9999908 , 0.91607307, 0.84471905, 0.7833553 , 0.73005756,
   0.68336001, 0.64212767, 0.60546881, 0.57267363, 0.5431705 ,
   0.51649874])