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.

Improve this question
11
  • 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. Commented Nov 10, 2022 at 15:25
  • 1
    If you give the code it makes easier to give you an answer. Commented Nov 11, 2022 at 19:29
  • 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. Commented Nov 14, 2022 at 21:16
  • 1
    Aren't you struck by the possible number of solutions ? Commented Nov 15, 2022 at 12:49
  • 1
    Did I say that it was impossible ? Commented Nov 15, 2022 at 15:14

1 Answer 1

Reset to default
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])
Improve this answer
Sign up to request clarification or add additional context in comments.

Comments

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.