Solving system of stiff ODE's using multivariate Newton Raphson method in Python

Aim:-> To solve ODE's using multivariant Newton Raphson method

Introduction:-

The Newton Raphson method is used to get the solutions of nonlinear equations system. The Newton-Raphson method (also known as Newton's method) is a way to find an approximation for the root of a real-valued function `f(x)=0`. The methods rely on the idea that a continuous and differentiable function can be approximated by a straight-line tangent to it.

Following ODE's are solved using Newton Raphson method

`(dy_1)/dt = -0.04y_1 + 10^4*y_2*y_3`

`(dy_2)/dt = 0.04y_1 -10^4*y_2*y_3 -3xx10^7*y_2^2`

`(dy_3)/dt = 3xx10^7*y_2^2`

As we have initial guess values for `y_1,y_2,y_3`

`y_1=1`

`y_2=0`

`y_3=0`

It is visible in above equation that the values of `y_1,y_2,y_3` may change rapidly due to nonlinear terms with higher order which are present in ODE's i.e. above equations shows stiff ODE's. So using implicit method is best suitable option for solving such stiff ODE's.

For all of the below equation superscript 'n' shows new node or new timestep while 'n-1' shows previous or old timestep or node.

By using FDM,

`(y^(n)-y^(n-1))/(Delta t)=C`

`y^(n)-(Delta txxC)=y^(n-1)`

So above ODE's can be calculate as,

`(y_1^n - y_1^(n-1))/(Delta t) = (-0.04y_1 + 10^4*y_2*y_3)^n`

`(y_2^n - y_2^(n-1))/(Delta t) = (0.04y_1 -10^4*y_2*y_3 -3xx10^7*y_2^2)^n`

`(y_3^n - y_3^(n-1))/(Delta t) = (3xx10^7*y_2^2)^n`

Now rearrange above equation.

Hence we get system of linear equations,

`F_1^n = y_1^(n-1) + (Delta t)(-0.04y_1 + 10^4*y_2*y_3)^n - y_1^n =0`

`F_2^n = y_2^(n-1) + (Delta t)(0.04y_1 -10^4*y_2*y_3 -3xx10^7*y_2^2)^n - y_2^n =0`

`F_3^n = y_3^(n-1) + (Delta t)(3xx10^7*y_2^2)^n - y_3^n =0 `

Numerically Jacobian can be calculated as-

`(df)/(dx)=(f(x+h)-f(x))/h`

Hence Jacobian matrix is calculated as-

`J=[[(df_1)/(dy_1),(df_1)/(dy_2),(df_1)/(dy_3)],[(df_2)/(dy_1),(df_2)/(dy_2),(df_2)/(dy_3)],[(df_3)/(dy_1),(df_3)/(dy_2),(df_3)/(dy_3)]]`

So now by using Newton Raphson method the value can be calculated as,

`Y_("new") = Y_("guess")-alphaJ^(-1)F`

Putting values to get final equation-

`y^("new")=[[y_1^("new")],[y_2^("new")],[y_3^("new")]]=[[y_1^("old")],[y_2^("old")],[y_3^("old")]]-alpha[[(df_1)/(dy_1),(df_1)/(dy_2),(df_1)/(dy_3)],[(df_2)/(dy_1),(df_2)/(dy_2),(df_2)/(dy_3)],[(df_3)/(dy_1),(df_3)/(dy_2),(df_3)/(dy_3)]]^(-1)[[F_1],[F_2],[F_3]]`

Python code:-

'''
Code for Multi varient Newton Raphson method

dy1/dt = -0.04*y1 + 10^4*y2*y3
dy2/dt = 0.04*y1 -10^4*y2*y3 -3*10^7*y2^2
dy3/dt = 3*10^7*y2^2 

initial conditions 
y1 = 1; y2 = 0; y3 = 0

relaxation factor = alpha = 1
'''

# Importing libraries
import numpy as np 
from numpy.linalg import inv
import matplotlib.pyplot as plt 

# Initial timing
t_init = 0 
# Total time
total_time = 10 # minutes
t_final = total_time*60 #seconds
dt = 0.5
# Number of nodes
n = int((t_final - t_init)/dt)
# Relaxation factor
alpha = 1

# Functions to solve equations

# For 1st equation
def f1(y1,y2,y3,dt,y1_old):
	return (y1_old + dt*(-0.04*y1 + y2*y3*pow(10,4)) - y1) 

# For 2nd equation
def f2(y1,y2,y3,dt,y2_old):
	return (y2_old + dt*(0.04*y1 - y2*y3*pow(10,4) - 3*pow(10,7)*pow(y2,2)) - y2)

# For 3rd equation
def f3(y1,y2,y3,dt,y3_old):
	return (y3_old + dt*(3*pow(10,7)*pow(y2,2)) - y3)


# Function to find Jacobian matrix for provided equations
def jacobi(y1,y1_old,y2,y2_old,y3,y3_old,dt):

	h = 1e-6

	# Jacobian matrix
	j = np.ones((3,3))

	#1st row
	j[0,0] = (f1(y1+h,y2,y3,dt,y1_old) - f1(y1,y2,y3,dt,y1_old))/h
	j[0,1] = (f1(y1,y2+h,y3,dt,y1_old) - f1(y1,y2,y3,dt,y1_old))/h
	j[0,2] = (f1(y1,y2,y3+h,dt,y1_old) - f1(y1,y2,y3,dt,y1_old))/h

	#2nd row
	j[1,0] = (f2(y1+h,y2,y3,dt,y2_old) - f2(y1,y2,y3,dt,y2_old))/h
	j[1,1] = (f2(y1,y2+h,y3,dt,y2_old) - f2(y1,y2,y3,dt,y2_old))/h
	j[1,2] = (f2(y1,y2,y3+h,dt,y2_old) - f2(y1,y2,y3,dt,y2_old))/h

	#3rd row
	j[2,0] = (f3(y1+h,y2,y3,dt,y3_old) - f3(y1,y2,y3,dt,y3_old))/h
	j[2,1] = (f3(y1,y2+h,y3,dt,y3_old) - f3(y1,y2,y3,dt,y3_old))/h
	j[2,2] = (f3(y1,y2,y3+h,dt,y3_old) - f3(y1,y2,y3,dt,y3_old))/h

	return j

# Initial condition at t_init
y1_old = 1 
y2_old = 0
y3_old = 0

y_old = np.ones((3,1))
y_old[0] = y1_old
y_old[1] = y2_old
y_old[2] = y3_old

F = np.copy(y_old)

# Tolerence
tol = 1e-6

Total_iter = 0

# Empty arrays creation
y_1 = []
y_2 = []
y_3 = []
it_mt = []

# Time loop
for i in range(0,n):

	iter = 1
	error = 9e9

	# Space loop
	while tol < error:

		# Calling Jacobi function
		j = jacobi(y_old[0],y1_old,y_old[1],y2_old,y_old[2],y3_old,dt)

		F[0] = f1(y_old[0],y_old[1],y_old[2],dt,y1_old)
		F[1] = f2(y_old[0],y_old[1],y_old[2],dt,y2_old)
		F[2] = f3(y_old[0],y_old[1],y_old[2],dt,y3_old)

		# Finit new guess values for y
		y_new = y_old - alpha*np.matmul(inv(j), F)

		# Calculating error
		error = np.max(np.abs(y_new - y_old))

		# Updating old values (Previous guess values)
		y_old = y_new
		
		# Updating iteration value
		iter = iter + 1

	# End of space loop

	# Printing of matrix after getting coverged values for Y
	log_msg = 'Itereation = {0} y1 = {1} y2 = {2} y3 = {3}'.format(iter, y1_old, y2_old, y3_old)
	print(log_msg)

	# Calculate total number of iterations
	Total_iter = Total_iter + iter

	# Updating old values (Previous time step values)
	y_old = y_new
	y1_old = y_new[0]
	y2_old = y_new[1]
	y3_old = y_new[2]

	y_1.append(y_new[0])
	y_2.append(y_new[1])
	y_3.append(y_new[2])
	it_mt.append(iter)

# End of time loop

# Printing total numbers of iteration
print('Total number of iterations are =',Total_iter)

# Plotting
plt.plot(y_1,'--',linewidth=1.5,color='purple')
plt.plot(y_2,'-.',linewidth=1.5,color='red')
plt.plot(y_3,linewidth=1.5,color='green')
plt.xlabel('Number of nodes')
plt.ylabel('Species')
plt.grid ('on')
plt.legend(['y_1', 'y_2', 'y_3'])
plt.title('For dt = '+str(dt)+' sec with number of '+str(Total_iter)+' iterations & '+str(alpha)+' relaxation factor')
plt.show()

Output plots and data:-

Output for final 10 time steps:- Itereation = 3 y1 = [0.40123628] y2 = [2.64538116e-06] y3 = [0.59876107] Itereation = 3 y1 = [0.4011314] y2 = [2.64424739e-06] y3 = [0.59886595] Itereation = 3 y1 = [0.40102661] y2 = [2.64311496e-06] y3 = [0.59897074] Itereation = 3 y1 = [0.40092191] y2 = [2.64198388e-06] y3 = [0.59907544] Itereation = 3 y1 = [0.4008173] y2 = [2.64085414e-06] y3 = [0.59918005] Itereation = 3 y1 = [0.40071278] y2 = [2.63972573e-06] y3 = [0.59928458] Itereation = 3 y1 = [0.40060835] y2 = [2.63859866e-06] y3 = [0.59938901] Itereation = 3 y1 = [0.40050401] y2 = [2.63747292e-06] y3 = [0.59949335] Itereation = 3 y1 = [0.40039975] y2 = [2.63634851e-06] y3 = [0.59959761] Itereation = 3 y1 = [0.40029559] y2 = [2.63522542e-06] y3 = [0.59970177] Total number of iterations are = 3669

Output for final 10 time steps:- Itereation = 3 y1 = [0.40057794] y2 = [2.63827047e-06] y3 = [0.59941943] Itereation = 3 y1 = [0.40052576] y2 = [2.6377075e-06] y3 = [0.59947161] Itereation = 3 y1 = [0.4004736] y2 = [2.63714487e-06] y3 = [0.59952377] Itereation = 3 y1 = [0.40042146] y2 = [2.63658257e-06] y3 = [0.5995759] Itereation = 3 y1 = [0.40036935] y2 = [2.63602061e-06] y3 = [0.59962802] Itereation = 3 y1 = [0.40031726] y2 = [2.63545897e-06] y3 = [0.59968011] Itereation = 3 y1 = [0.40026519] y2 = [2.63489767e-06] y3 = [0.59973218] Itereation = 3 y1 = [0.40021314] y2 = [2.63433669e-06] y3 = [0.59978423] Itereation = 3 y1 = [0.40016111] y2 = [2.63377605e-06] y3 = [0.59983625] Itereation = 3 y1 = [0.40010911] y2 = [2.63321574e-06] y3 = [0.59988826] Total number of iterations are = 7250

Output for final 10 time steps:- Itereation = 3 y1 = [0.40255859] y2 = [2.65970782e-06] y3 = [0.59743875] Itereation = 3 y1 = [0.40234674] y2 = [2.65740848e-06] y3 = [0.5976506] Itereation = 3 y1 = [0.40213526] y2 = [2.65511464e-06] y3 = [0.59786209] Itereation = 3 y1 = [0.40192413] y2 = [2.65282628e-06] y3 = [0.59807321] Itereation = 3 y1 = [0.40171338] y2 = [2.65054337e-06] y3 = [0.59828397] Itereation = 3 y1 = [0.40150298] y2 = [2.6482659e-06] y3 = [0.59849437] Itereation = 3 y1 = [0.40129294] y2 = [2.64599383e-06] y3 = [0.59870441] Itereation = 3 y1 = [0.40108327] y2 = [2.64372716e-06] y3 = [0.59891409] Itereation = 3 y1 = [0.40087395] y2 = [2.64146585e-06] y3 = [0.59912341] Itereation = 3 y1 = [0.40066499] y2 = [2.63920989e-06] y3 = [0.59933237] Total number of iterations are = 1893

Output for final 10 time steps:- Itereation = 6 y1 = [0.40120493] y2 = [2.64504218e-06] y3 = [0.59879243] Itereation = 6 y1 = [0.40110004] y2 = [2.64390845e-06] y3 = [0.59889731] Itereation = 6 y1 = [0.40099525] y2 = [2.64277606e-06] y3 = [0.59900211] Itereation = 6 y1 = [0.40089054] y2 = [2.64164502e-06] y3 = [0.59910682] Itereation = 6 y1 = [0.40078592] y2 = [2.64051532e-06] y3 = [0.59921144] Itereation = 5 y1 = [0.40068139] y2 = [2.63938695e-06] y3 = [0.59931597] Itereation = 5 y1 = [0.40057716] y2 = [2.63826183e-06] y3 = [0.59942021] Itereation = 5 y1 = [0.40047301] y2 = [2.63713828e-06] y3 = [0.59952436] Itereation = 5 y1 = [0.40036894] y2 = [2.63601605e-06] y3 = [0.59962842] Itereation = 5 y1 = [0.40026497] y2 = [2.63489515e-06] y3 = [0.59973239] Total number of iterations are = 7452

Output for final 10 time steps:- Itereation = 5 y1 = [0.40132481] y2 = [2.64633853e-06] y3 = [0.59867255] Itereation = 5 y1 = [0.40122002] y2 = [2.64520544e-06] y3 = [0.59877733] Itereation = 5 y1 = [0.40111532] y2 = [2.6440737e-06] y3 = [0.59888203] Itereation = 5 y1 = [0.40101072] y2 = [2.64294329e-06] y3 = [0.59898664] Itereation = 5 y1 = [0.4009062] y2 = [2.64181423e-06] y3 = [0.59909116] Itereation = 5 y1 = [0.40080177] y2 = [2.6406865e-06] y3 = [0.59919559] Itereation = 5 y1 = [0.40069743] y2 = [2.6395601e-06] y3 = [0.59929993] Itereation = 5 y1 = [0.40059318] y2 = [2.63843503e-06] y3 = [0.59940419] Itereation = 5 y1 = [0.40048901] y2 = [2.63731128e-06] y3 = [0.59950835] Itereation = 5 y1 = [0.40038494] y2 = [2.63618886e-06] y3 = [0.59961242] Total number of iterations are = 6971

Conclusion:- From all above output plots and data it is observed and can be conclude that-

1. As we move forward in time the value of functions varies rapidly.
2. Change in time step doesn't make any big difference in final solution.
3. Small value of time step takes longer time (iterations) to converged result which may increase cost.
4. Considering small timestep increases number of nodes which gives more precise solution curve (in plot).
5. As implicit scheme is considered we get stable solution but setting implicit scheme is more complex and hectic.
6. For over relaxed setting number of iterations are higher than under relaxed setting.