Solve ODE using Multivariate Newton Raphson Method

Objective: Solving a set of ODE equations using the multivariate Newton Raphson method

Below are the given equations

As per stiff problem, it is better to use implicit backward differentiation because the problem solution doesn't get blown although the accuracy can affected 

J is J-1⋅F</math>" role="presentation">J−1⋅F<math xmlns="http://www.w3.org/1998/Math/MathML">J-1⋅F</math>[removed]J^-1*F[removed] jacoboan matrix with size of function rows and variable columns

F is the function i take at the start of the problem

f1=( -0.04*y1+ 1e4*y2*y3 ) * dt - y1

f2=y2_old+(0.04*y1-1e4*y2*y3-3*1e7*y2**2)*dt - y2

f3=y3_old+(3*1e7*y2**2)*dt-y3

Here, jacobian is calculated numerically by using the Newton-Raphson method

#As per the below code

Matplotlib, numpy, and numpy.linalg (to perform inverse calculation for Jacobi matrix) are imported

As i engage in a 10-minute simulation, one of the pivotal steps is to meticulously select a stable time step. The choice of time step is of paramount importance, particularly in the case of an explicit method. However, in our current scenario, an implicit approach is employed, thus rendering the specific time step less critical as the solution will not destabilize. Nonetheless, it's essential to note that the accuracy of the outcomes may be impacted. Therefore, it is advisable to opt for a smaller time step to ensure precision in results.

Here, assume initial conditions for

Y1=1

Y2=0

Y3=0

and alpha=1

 # Solve multivarite newton Raphson method

import matplotlib.pyplot as plt
import numpy as np
from numpy.linalg import inv

t_start=0
t_stop=600 
dt=0.01

def f1(y1_old,y1,y2,y3,dt):
 return y1_old + ( -0.04*y1+ 1e4*y2*y3 ) * dt - y1

def f2(y2_old,y1,y2,y3,dt):
 return y2_old+(0.04*y1-1e4*y2*y3-3*1e7*y2**2)*dt - y2

def f3(y3_old,y1,y2,y3,dt):
 return y3_old+(3*1e7*y2**2)*dt-y3

def jac(Y_old,Y_guess,dy,dt):
 
 J=np.ones((3,3))

 J[0,0]= (f1(Y_old[0],Y_guess[0]+dy,Y_guess[1],Y_guess[2],dt)-f1(Y_old[0],Y_guess[0],Y_guess[1],Y_guess[2],dt))/dy
 J[0,1]= (f1(Y_old[0],Y_guess[0],Y_guess[1]+dy,Y_guess[2],dt)-f1(Y_old[0],Y_guess[0],Y_guess[1],Y_guess[2],dt))/dy
 J[0,2]= (f1(Y_old[0],Y_guess[0],Y_guess[1],Y_guess[2]+dy,dt)-f1(Y_old[0],Y_guess[0],Y_guess[1],Y_guess[2],dt))/dy
 J[1,0]= (f2(Y_old[1],Y_guess[0]+dy,Y_guess[1],Y_guess[2],dt)-f2(Y_old[1],Y_guess[0],Y_guess[1],Y_guess[2],dt))/dy
 J[1,1]= (f2(Y_old[1],Y_guess[0],Y_guess[1]+dy,Y_guess[2],dt)-f2(Y_old[1],Y_guess[0],Y_guess[1],Y_guess[2],dt))/dy
 J[1,2]= (f2(Y_old[1],Y_guess[0],Y_guess[1],Y_guess[2]+dy,dt)-f2(Y_old[1],Y_guess[0],Y_guess[1],Y_guess[2],dt))/dy
 J[2,0]= (f3(Y_old[2],Y_guess[0]+dy,Y_guess[1],Y_guess[2],dt)-f3(Y_old[2],Y_guess[0],Y_guess[1],Y_guess[2],dt))/dy
 J[2,1]= (f3(Y_old[2],Y_guess[0],Y_guess[1]+dy,Y_guess[2],dt)-f3(Y_old[2],Y_guess[0],Y_guess[1],Y_guess[2],dt))/dy
 J[2,2]= (f3(Y_old[2],Y_guess[0],Y_guess[1],Y_guess[2]+dy,dt)-f3(Y_old[2],Y_guess[0],Y_guess[1],Y_guess[2],dt))/dy
 
 return J

t=np.arange(t_start,t_stop,dt)
T=len(t)
y1=1
y2=0
y3=0
Y_old=np.ones((3,1)) 
Y_guess=np.ones((3,1))
F=np.copy(Y_old)
dy=1e-5

iter = 0
Y_old[0]=y1
Y_old[1]=y2
Y_old[2]=y3


Y_guess[0]=y1
Y_guess[1]=y2
Y_guess[2]=y3
Y1_Final=[]
Y2_Final=[]
Y3_Final=[]

for i in range(0,T):
  
  
  error=9e9
  alpha=1
  tol=1e-3
  while error > tol:

   J=jac(Y_old,Y_guess,dy,dt);

   F[0]=f1(Y_old[0],Y_guess[0],Y_guess[1],Y_guess[2],dt);
   F[1]=f2(Y_old[1],Y_guess[0],Y_guess[1],Y_guess[2],dt);
   F[2]=f3(Y_old[2],Y_guess[0],Y_guess[1],Y_guess[2],dt);

   Y_new = Y_guess - alpha*np.matmul(inv(J),F)
   error=np.max(np.abs(Y_new-Y_guess))
   Y_old=np.copy(Y_new)
   Y_guess=np.copy(Y_new)
   
 

  iter = iter + 1
  log_message=\'iteration # = {0} y1={1} y2={2} y3={3} error={1}\'. format(iter,Y_new[0],Y_new[1],Y_new[2],error)
  print(log_message)
  Y1_Final.append(Y_new[0]) 
  Y2_Final.append(Y_new[1])
  Y3_Final.append(Y_new[2]) 


plt.plot(t,Y1_Final)
plt.plot(t,Y2_Final)
plt.plot(t,Y3_Final)
plt.legend([\'Y1\',\'Y2\',\'Y3\'])
plt.xlabel(\'Time(s)\')
plt.ylabel(\'Y\')

plt.show()

After completing the simulation for the time i assumed let us plot a graph(Time vs variable(Y1,Y2,Y3)) to see the profiles. It looks like the solution is fine

Here, the Y axis is plotted with the Y1, Y2, and Y3 with respect to time