Week 6 - Multivariate Newton Rhapson Solver

In the project, we are going to learn the idea of how the rate of change of production of different species which are in ODE  form for each reaction. These ODE's are set of coupled equation ie the rate of one specie is dependent on the on its concentration as well on others. In order to solve these ODE, we use one of the numerical method called Multivariate Newton Raphston methods. The Multivariate Newton Raphston method finds the roots of the function iteratively.

`Y_n =Y_o - (Inv(J) f)`

Here instead of species, we will use coupled variables 

`(dy_1)/(dt)=-0.04y_1+10^4y_2y_4`

`dy_2/dt = 0.04y_1-10^4y2y_3 -3.10^7y_2`

`dy_3/dt=3.10^7y_2^2`

The above are the given ODE's

The RHS side which is time differential is discretized using backward differencing scheme along with the implicit method is implemented and then further the system of ODE's are rearranged that is according to Multivariate Newton Raphston method to find the roots

Discretization of Time differentials

`(dy_1^n-dy_1^(n-1))/(dt)=-0.04y_1+10^4y_2y_4`

`(dy_2^n-dy_2^(n-1))/dt = 0.04y_1-10^4y2y_3 -3.10^7y_2^2`

`(dy_3^n-dy_3^(n-1))/dt=3.10^7y_2^2`

Rearrangement

`f_1=y_1^n+h(0.04(y_1)^n-10^4(y_2)^n(y_4)^n)-y_1^(n-1)`

`f_2=y_2^n -h(0.04(y_1)^n-10^4(y2)^n(y_3)^n -3.10^7(y_2^2)^n)-y_2^(n-1) =0`

`y_3^n=h(3.10^7(y_2^2)^n)-y_3^(n-1)=0`

The jacobian form 

The differentials in the jacobians are calculated using the forward differencing

`(df_i)/dy_i=(f_i(y_i+triangley_i))/(triangley_i)`

and here `triangley_i` is taken as small  and the value is taken is1e-8

The Python Code 

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


#Defining the root solving function
def f1(y1_n,y2_n,y3_n,y1_o,h):
	return y1_n + 0.04*h*y1_n - 1e4*h*y2_n*y3_n - y1_o

def f2(y1_n,y2_n,y3_n,y2_o,h):
	return y2_n -0.04*h*y1_n + 1e4*h*y2_n*y3_n + 3e7*h*pow(y2_n,2) - y2_o

def f3(y1_n,y2_n,y3_n,y3_o,h):
	return y3_n - 3e7*h*pow(y2_n,2) - y3_o

#Calculating Jacobian using spcial forward differencing
def jac(Y_n,Y_o,h):
	dy = 1e-8
	y1_n = Y_n[0]
	y2_n = Y_n[1]
	y3_n = Y_n[2]

	y1_o = Y_o[0]
	y2_o = Y_o[1]
	y3_o = Y_o[2]

	#Jacobian matrix for all the three equation and three unkown
	J = np.zeros((3,3))

	J[0,0] = (f1(y1_n+dy,y2_n,y3_n,y1_o,h) - f1(y1_n,y2_n,y3_n,y1_o,h))/dy
	J[0,1] = (f1(y1_n,y2_n+dy,y3_n,y1_o,h) - f1(y1_n,y2_n,y3_n,y1_o,h))/dy
	J[0,2] = (f1(y1_n,y2_n,y3_n+dy,y1_o,h) - f1(y1_n,y2_n,y3_n,y1_o,h))/dy

	J[1,0] = (f2(y1_n+dy,y2_n,y3_n,y2_o,h) - f2(y1_n,y2_n,y3_n,y2_o,h))/dy
	J[1,1] = (f2(y1_n,y2_n+dy,y3_n,y2_o,h) - f2(y1_n,y2_n,y3_n,y2_o,h))/dy
	J[1,2] = (f2(y1_n,y2_n,y3_n+dy,y2_o,h) - f2(y1_n,y2_n,y3_n,y2_o,h))/dy

	J[2,0] = (f3(y1_n+dy,y2_n,y3_n,y3_o,h) - f3(y1_n,y2_n,y3_n,y3_o,h))/dy
	J[2,1] = (f3(y1_n,y2_n+dy,y3_n,y3_o,h) - f3(y1_n,y2_n,y3_n,y3_o,h))/dy
	J[2,2] = (f3(y1_n,y2_n,y3_n+dy,y3_o,h) - f3(y1_n,y2_n,y3_n,y3_o,h))/dy

	return J

#Defining the initial condition
Y_o = np.zeros((3,1))
Y_o[0] = 1

Y_n = np.zeros((3,1))

F= np.copy(Y_o)

error = 1e9
alpha = 1
tol = 1e-10
iter = 0
t= np.arange(0,600,step=0.1)
h=t[1] - t[0]

y1 = [1]*len(t)
y2 = [0]*len(t)
y3 = [0]*len(t)
iteration=[]

# Outer time loop
for k in range(1,len(t)):
	y1_o = Y_o[0]
	y2_o = Y_o[1]
	y3_o = Y_o[2]
	Y_fm  = Y_o

	#Inner iteratve solver loop
	while error >= tol:
		J = jac(Y_n,Y_o,h)
		y1_fm = Y_fm[0]
		y2_fm = Y_fm[1]
		y3_fm = Y_fm[2]

		F[0] = f1(y1_fm,y2_fm,y3_fm,y1_o,h)
		F[1] = f2(y1_fm,y2_fm,y3_fm,y2_o,h)
		F[2] = f3(y1_fm,y2_fm,y3_fm,y3_o,h)

		Y_n = Y_fm - alpha*np.matmul(inv (J), F)
		error = max(abs(Y_n-Y_fm))

		Y_fm = Y_n



	Y_o = Y_n
	y1[k] = Y_n[0]
	y2[k] = Y_n[1]
	y3[k] = Y_n[2]

	print(f'time ={round(t[k], 1)} sec, iteration No = {k},y1 = {round(Y_n[0][0], 2)}, y2 = {round(Y_n[1][0], 2)}, y3 = {round(Y_n[2][0], 2)}')


	iter =iter+1
	iteration.append(iter)
	error=1e9
plt.plot(t,y1,color = 'r')
plt.plot(t,y2,color = 'b')
plt.plot(t,y3,color = 'g')
plt.xlabel('Time(sec)')
plt.ylabel('Variable')
plt.legend(['y1','y2','y3'])
plt.show()

Results

Time Step=0.1

Time Step=0.001

Time step =1

Since the implicit scheme is being used, the solution will stable regardless of the time step being used.but what needs to consider to choose an ideal time step is an accuracy and No of iteration.The time step of  0.01 took too many iterations which takes a large amount of memory space while the step of 1 took less number of iteration about 599 and we cannot conclude its accuracy. The step time step of 0.1 shown results somewhere between the two, having a moderate amount of iteration and not taking much memory space