1

enter image description here

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.

simon
  • 83
  • 1
  • 9
  • 1
    Yes, you can use iteration using the old values inside the matrix. Or you treat the whole as non-linear system and apply Newton's method. This will give a slightly different fixed-point iteration that hopefully converges faster. – Lutz Lehmann Nov 10 '22 at 15:25
  • 1
    If you give the code it makes easier to give you an answer. – Bob Nov 11 '22 at 19:29
  • Thanks! I will give it a go (Lutz Lehmann). The reason why I added zero code, is because I had no idea what maths to use to tackle the question so any code would have been useless. But now I know its iterative method that will help I will now be able to, once done I will post the solution on here (Bob). – simon Nov 12 '22 at 11:07
  • 1
    Notice that your system is quadratic, with a special structure. It reads `p² + T p = R`, where `p²` is a vector made of the squares of the components of `p` and `T` is tridiagonal. It potentially has 2^11 solutions. –  Nov 14 '22 at 21:16
  • would newtons method be able to solve this or would another method be needed @YvesDaoust? – simon Nov 15 '22 at 12:20
  • 1
    Aren't you struck by the possible number of solutions ? –  Nov 15 '22 at 12:49
  • Of course but if there someone else can do it this means it must be possible @YvesDaoust. Ill keep trying feel like Im getting closer with more information – simon Nov 15 '22 at 14:08
  • 1
    Did I say that it was impossible ? –  Nov 15 '22 at 15:14
  • Never said you said it was impossible. But your previous comment, conveyed that this is a difficult problem, I'm just saying its difficult but possible. @YvesDaoust – simon Nov 15 '22 at 16:06
  • 1
    Do you just don't care that there can be many solutions ? –  Nov 15 '22 at 16:14
  • Yep, @YvesDaoust – simon Nov 15 '22 at 18:30

1 Answers1

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])
simon
  • 83
  • 1
  • 9