Implementing multivariate Newton-Raphson method

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_text{new} = y_text{guess} - alpha {f(y_text{guess})}/{f'(y_text{guess})}`

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:

`bb{Y_text{new}} = bb{Y_text{guess}} - alpha bb{J^{-1}F}`

where `bb{J}` is the Jacobian matrix of the system of coupled ODEs and `alpha` is the relaxation parameter.

FORMULATING THE SYSTEM:

For this project, we will be solving a system of coupled ODEs given as:

`{dy_1}/{dt} = -0.04 y_1 + 10^4 y_2 y_3`

`{dy_2}/{dt} = 0.04 y_1 - 10^4 y_2 y_3 - 3 times 10^7 y_2^2`

`{dy_3}/{dt} = 3times10^7 y_2^2`

with the initial conditions as: `y_1(0) = 1, y_2(0) = 0` and `y_3(0) = 0`

The time-differentials are discretized using the backward difference scheme and implicit method is implemented for a stable solution.

`{y_1^n - y_1^{n-1}}/{Delta t} = -0.04 y_1^n + 10^4 y_2^n y_3^n`

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

`{y_3^n - y_3^{n-1}}/{Delta t} = 3 times 10^7 (y_2^n)^2`

The superscript `n` denotes the time-step at which the values are taken. To solve the system using the multivariate Newton-Raphson solver, the system needs to be converted to a system of root-finding functions. Upon rearranging the terms in the above equations and taking `h = Delta t`, we get:

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

`f_2^n = y_2^n - h(0.04 y_1^n - 10^4 y_2^n y_3^n - 3 times 10^7 (y_2^n)^2) - y_2^{n-1}=0`

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

Now that we have the root finding functions, we can calculate the Jacobian matrix as:

`bb{J} = [({partial f_1^n}/{partial y_1^n},,{partial f_1^n}/{partial y_2^n},,{partial f_1^n}/{partial y_3^n}),({partial f_2^n}/{partial y_1^n},,{partial f_2^n}/{partial y_2^n},,{partial f_2^n}/{partial y_3^n}),({partial f_3^n}/{partial y_1^n},,{partial f_3^n}/{partial y_2^n},,{partial f_3^n}/{partial y_3^n})]`

The system can then be solved for `y_1, y_2` and `y_3` as:

`bb{Y_{k+1}^n} = bb{Y_k^n} - alpha [bb{J^{-1} F}]^n`

Here, the subscript `k` indicates the iteration number within the time-step loop and `bb{F}` is the function matrix `bb{F} = [(f_1),(f_2),(f_3)]`.

IMPLEMENTATION IN PYTHON:

The solver is implemented using Python. The Jacobian matrix is calculated using numerical approximation instead of analytical differentiation. The differentials are calculated using forward differencing as:

`{partial f_i}/{partial y_i} = {f_i(y_i + Delta y_i) - f_i(y_i)}/{Delta y_i}`

Here, `Delta y_i` is taken to be as small as possible for accuracy: `Delta y_i = 10^{-8}` for `i = {1, 2, 3}`

The system of ODEs is solved for a physical time of `10 min` `(600` `sec)`. Since we are implementing the implicit Euler method, the solution will be stable irrespective of the time-step chosen. However, the accuracy of the solution depends highly on the time-step and hence the time-step is chosen to be `0.1` `sec`. The tolerance for the iterative error is set to `10^{-12}`. The relaxation parameter is taken to be `alpha = 1`. The Python code and the results obtained are shown below:

"""
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: Shrey Shah
Date: 30-Sep-19
"""

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='red',linewidth=2,label='$y_1$')
plt.plot(t,y2,color='black',linewidth=2,label='$y_2$')
plt.plot(t,y3,color='blue',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='black')
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()