Week 6 - Multivariate Newton Rhapson Solver

Aim: Code a multivariate Newton-Raphson Iterator

Objective:

• Solve the given set of differential equation by coding a multivariate newton-raphson iterator
• Solve using Implicit Euler method / Backward differencing for a given initial condition.

Theory:

In combustion mechanism, the rates of reactant comsumption and product formation are calculated in the form of ODEs. Hence, DEs are formed for each individual species. The set of these DEs are called as a system of coupled DE. These cannot be solved directly through integration, hence numerical methods are used to integrate and solve the system. One such method is the multivariate Newton-Raphson method, which is extended part of Newton-Raphson. The NR method solves for the root of a function through iterative technique (using guess values) as shown:

`y_{n e w} = y_{guess} - alpha((f(y_{guess}))/(f'(y_{guess})))`

The same concept when applied to a multivariate system results in:

`Y_{n e w} = Y_{guess} - alphabb{J}^-1bb{F}`

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

Deriving the system:

The given system of coupled ODEs are:

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

`(dy_2)/(dt) = 0.04y_1 + 10^4y_2y_3 - 3 times 10^7y_2^2`

`(dy_2)/(dt) = 3 times 10^7y_2^2`

with initial conditions, `y_1(0) = 1, y_2(0) = 0 a n d 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})/(Deltat) = -0.04y_1^n + 10^4y_2^ny_3^n`

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

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

where, `n` denotes the time-step at which the values are taken.To solve the system using NR solver, the system needs to be arranged in the form of root finding functions. After rearranging and replacing `Deltat` by `h`, we have,

`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 times 10^7(y_2^n)^2) - y_1^{n-1} = 0`

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

And the Jacobian is:

`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 be solved for `y_1`, `y_2` and `y_3` as:

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

where, `k` is the iteration number and `bb{F}` is the function matrix

`bb{F} = [[f_1], [f_2], [f_3]]`

Code:

The Python code for the NR solver is as shown. The Jacobian matrix is caculated by numerical approximation (forward differencing) as:

`(partial f_i)/(partial y_i) = (f_i(y_i + Deltay_i) - f_i(y_i))/(Deltay_i)`

where, `Deltay_i` is taken to be `10^-8` (smaller value gives greater accuracy)

The system of ODEs is solved for a time of 10min (600 sec). And since we have implemented implicit Euler method, the solution will be stable irrespctive of time-chosen. BUt, to obtain the solution, we need to provide a time step, and a time step of 0.1 sec was chosen. The tolerance for the iteration was set at `10^-12` and relaxation factor was set to 1.

"""
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
"""

# Importing Required Libraries
from numpy import *
from matplotlib.pyplot import *
from numpy.linalg import inv

# 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 = 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 = zeros((3,1))
Y_o[0] = 1

Y_n = zeros((3,1))
F = copy(Y_o)

# Defining the Newton-Raphson solver parameters
err = 1e9
alpha = 1
tol = 1e-12
count = 0
t = arange(0,600.00001,step=0.01)
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*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}\n'.format(round(t[k],1),round(Y_n[0][0],3),round(Y_n[1][0],3),round(Y_n[2][0],3))
	file = open("log_message.txt","a")
	file.write(log_message)
	# Resetting the error criteria
	count = 1
	err = 1e9

# Plotting the variables at each time-step
fig1, ax1 = subplots()
plot(t,y1,color='red',linewidth=2,label='$y_1$')
plot(t,y2,color='black',linewidth=2,label='$y_2$')
plot(t,y3,color='blue',linewidth=2,label='$y_3$')
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)
xlabel('Time (sec)',fontsize=14)
ylabel('$y_i$',fontsize=15)
title('Multivariate Newton-Raphson Solution for dt = '+str(h),fontsize=14,fontweight='bold')
legend(fontsize=16)
savefig("Solution.png")
show()

# Plotting the converged error for each time-step
fig2, ax2 = subplots()
semilogy(t,e,color='black')
grid('both',linestyle='--',linewidth=1)
xlabel('Time (sec)',fontsize=14)
ylabel('Max. Absolute Error',fontsize=15)
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)
savefig("Absolute_error.png")
show()
file.close()

Output:

Hence, solution of the coupled system of ODEs was obtained and the results were plotted and appended to a text file.

References:

• Computational Combustion using Python and Cantera - course provided by Skill Lync
• Anaconda - Anaconda Inc.
• JupyterLab - Project Jupyter
• Supporting Material: Higher Engineering Mathematics - BS Grewal