MULTIVARIATE NEWTON RHAPSON SOLVER FOR A SYSTEM OF COUPLED ODE s

INTRODUCTION

• While solving a combustion mechanism, the rates of production/depletion of different species are calculated in the form of ordinary differential equations (ODEs).

• The rate of production/dissipation of a particular species depends on the concentration of the species itself as well as the concentration of other species. Such set of ODEs is called a system of coupled differential equations.

• They cannot be solved directly through integration; hence numerical methods are used to integrate and solve the system of coupled differential equations. One such method is the multivariate Newton-Raphson method, which is an extension of the univariate Newton-Raphson method.

• The Newton-Raphson method solves for the root of a function (called the root-finding function) through iterative techniques. For a univariate system, we can find the root of the function through guess values as:

`y_n=y_(guess)-α*[f(y_(guess) )/(f^\' (y_(guess ) )]]`

Where `y_n` corrosponds to the new updated iteration value.

The above equation is iterated to the point where the difference between the new value and

the guess value reduces below a certain specified tolerance. This same method is

implemented to a multivariate system as:

`Y_n = Y_(guess) - alpha*(J^-1)* F`

where `J` is the Jacobian matrix of the system of coupled ODEs and `alpha` is the

relaxation parameter.

PROBLEM TO SOLVE

Initial Conditions: `y_1(0)= 1; y_2(0)=y_3(0)=0`

Now,

`(y_(1 )^n-y_(1 )^(n-1))/(∆t)=-0.04y_1^n+10^4y_2^ny_3^n`

`(y_(2 )^n-y_(2 )^(n-1))/(∆t)=0.04y_1^n-10^4y_2^ny_3^n -3*10^7(y_2^n)^2`

`(y_(3 )^n-y_(3 )^(n-1))/(∆t)=3*10^7(y_2^n)^2`

here, \'n\' is the time step at which values are taken.

Rearranging the terms in the above equation as a system of root finding functions, we get:

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

`f_2^n = y_2^n-h(0.04y_1^n-10^4y_2^ny_3^n-3*10^7(y_2^n)^2)-y_2^(n-1)=0`

`f_3^n = y_3^n-h(3*10^7(y_2^n)^2)-y_3^(n-1)=0`

In the above equations, `h=∆t`

We can calculate the Jacobian matrix for the above root finding equations as:

`J = [[(∂f_1^n)/(∂y_1^n),(∂f_1^n)/(∂y_2^n),(∂f_1^n)/(∂y_3^n)],[(∂f_2^n)/(∂y_1^n),(∂f_2^n)/(∂y_2^n),(∂f_2^n)/(∂y_3^n)],[(∂f_3^n)/(∂y_1^n),(∂f_3^n)/(∂y_2^n),(∂f_3^n)/(∂y_3^n)]]`

The system can now be solved for y1, y2 and y3 as:

`Y_(k+1)^n = Y_(k) ^n- alpha*[(J^-1)* F]^n`

k = Iteration number within the time step loop

F = function matrix. `F = [[f1],[f2],[f3]]`

PYTHON CODE

• `(∂f_i)/(∂y_i)=(f_i(y_i+∆y_i)-f_i(y_i))/(∆y_i)`
• `∆y_i=10^-8` for `i= {1,2,3}`
• Time step = 0.1 second
• Tolerance = `10^-12`
• `alpha=1`
• Total Simulation time = 600 s.

\"\"\"
Program to solve system of ODEs through multivariate Newton-Raphson method
System of ODEs:
y1\' = -0.04*y1 + 10^4*y2*y3
y2\' = 0.04*y1 - 10^4*y2*y3 - 3*10^7*y2*y2
y3\' = 3*10^7*y2*y2
Program by: Shouvik Bandopadhyay
\"\"\"

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

# Defining the root-solving functions
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 numerical Jacobian through spatial 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 3 equations and 3 unknowns
	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 conditions
Y_o = np.zeros((3,1))
Y_o[0] = 1
Y_n = np.zeros((3,1))
F = np.copy(Y_o)

# Defining the Newton-Raphson solver parameters
err = 1e9
alpha = 1
tol = 1e-12
count = 0
t = np.arange(0,600.00001,step=0.1)
h = t[1] - t[0]
y1 = [1]*len(t)
y2 = [0]*len(t)
y3 = [0]*len(t)
e = [1e-6]*len(t)

# 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_g = Y_o
	
	# Inner iterative solver loop
	while err >= tol:
		J = jac(Y_n,Y_o,h)
		y1_g = Y_g[0]
		y2_g = Y_g[1]
		y3_g = Y_g[2]
		F[0] = f1(y1_g,y2_g,y3_g,y1_o,h)
		F[1] = f2(y1_g,y2_g,y3_g,y2_o,h)
		F[2] = f3(y1_g,y2_g,y3_g,y3_o,h)
		Y_n = Y_g - alpha*np.matmul(inv(J),F)
		err = max(abs(Y_n - Y_g))
	
		# Updating the guess values for a new iteration
		Y_g = Y_n
		count += 1

	# Updating the new time-step values
	e[k] = err
	Y_o = Y_n
	y1[k] = Y_n[0]
	y2[k] = Y_n[1]
	y3[k] = Y_n[2]
	
	# Prompt message for each time-step
	log_message = \'Time = {0} sec, y1 = {1}, y2 = {2}, y3 = {3}\'.format
	(round(t[k],1),round(Y_n[0][0],3),round(Y_n[1][0],3),round(Y_n[2][0],3))
	print(log_message)
	
	# Resetting the error criteria
	count = 1
	err = 1e9

# Plotting the variables at each time-step
fig1, ax1 = plt.subplots()
plt.plot(t,y1,color=\'green\',linewidth=2,label=\'$y_1$\')
plt.plot(t,y2,color=\'purple\',linewidth=2,label=\'$y_2$\')
plt.plot(t,y3,color=\'brown\',linewidth=2,label=\'$y_3$\')
plt.grid(\'both\',linestyle=\'--\',linewidth=1)
xticks = [-100, 0, 100, 200, 300, 400, 500, 600, 700]
ax1.set_xticklabels(xticks,rotation=0,fontsize=12)
yticks = [-0.2, 0, 0.2, 0.4, 0.6, 0.8, 1, 1.2]
ax1.set_yticklabels(yticks,rotation=0,fontsize=12)
plt.xlabel(\'Time (sec)\',fontsize=14)
plt.ylabel(\'$y_i$\',fontsize=15)
plt.title(\'Multivariate Newton-Raphson Solution for dt = \'+str(h),fontsize=14,fontweight=\'bold\')
plt.legend(fontsize=16)
plt.show()

# Plotting the converged error for each time-step
fig2, ax2 = plt.subplots()
plt.semilogy(t,e,color=\'red\')
plt.grid(\'both\',linestyle=\'--\',linewidth=1)
plt.xlabel(\'Time (sec)\',fontsize=14)
plt.ylabel(\'Max. Absolute Error\',fontsize=15)
plt.title(\'Absolute error during each time-step for dt = \'+str(h),fontsize=14,fontweight=\'bold\')
xticks = [-100, 0, 100, 200, 300, 400, 500, 600, 700]
ax2.set_xticklabels(xticks,rotation=0,fontsize=12)
plt.show()

SIMULATION OUTPUTS

Result

y1 = 0.4; y2 = 0; y3 = 0.6