Multivariate Newton Rhapson Solver

INTRODUCTION :-

In order to estimate the combustion process inside a combustion chamber it is fundamental to evaluate the molar fractions of the gases which constitute the burned mixture. A non linear set of equations must be solved to evaluate the mixture composition.

This challenge deals with the study of solving the coupled ordinary differential equations involving about the thermodynamics of combustion, using a new approach to solve the set of non-linear equations system we obtain at each time step by Newton-Raphson method.

NEWTON RAPSHON APPROACH FOR MULTIVARIABLE :-

We have already solved for single variables which is written below mathematically in the below form :-

` x^(i+1)=x^i - α (f(x^i)/(J(x^i)))`where `J(x^i) = f\'(x^i)`

For a multivariable case all of the variables are converted to vector form and as we can divide it by J scalar quantity it now becomes and inverse jacobian matrix.

`X^(i+1)=X^i - α (J^-1F)` where J is jacobian matrix, F is Function matrix and α is the relaxation factor.

1.Solving the system :-

The system of coupled ODE\'s are given as :-

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

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

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

As this is coupled O.D.E i.e. yields a stiff system we use Implicit/Back-ward difference solving technique and the system becomes :-

`((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_2^n - y_2^(n-1))/(Δt) = 3.10^7(y_2^n)^2``

Substituting values and taking time step as h i.e. Δt

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

`y_2^n = y_2^(n-1) + h(0.04y_1^n - 10^4y_2^ny_3^n - 3.10^7(y_2^n)^2)`

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

Doing R.H.S-L.H.S yield us values with functions

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

`f_2= y_2^(n-1) + h(0.04y_1^n - 10^4y_2^ny_3^n - 3.10^7(y_2^n)^2)-y_2^n`

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

2. The function matrix can be written as follows :-

`F = ((f_1),(f_2),(f_3))`

3. The function matrix can be written as follows :-

As discussed above that Jacobian matrix is just the differential of function \'f\' at different variables available. Hence, for the given set of 3 equations and 3 variables. The jacobian matrix is 3x3 matrix.

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

4. Code implementation in python :-

The differentials that are to be solved are simplified as below with the general method of calculus differentitation given as :-

`(∂f)/(∂y) = (f(y+Δy) - f(y))/(Δy)`

Here Δy is taken as 1e-8. As we are implementing implicit method, solution shall be irresepective of time-step chosen as this method is quite stable. Regardless of that accuracy is dependent on small values of step size and we should chose them carefully and tolerance is taken to value of 1e-12 with relaxation factor as 1.

The python code is written below :-

\"Program to solve the system of non-linear coupled equations using Newton-Raphson method\"
\"\"\"
Equations are :-

1. dy1/dt = -0.04y1 + 10(power,4)*y2*y3
2. dy2/dt =  0.04y1 - 10(power,4)*y2*y3 -3*10(power,7)*y2(power,2)
3. dy3/dt =  3*10(power,7)*y2(power,2).

Modules required are :-

1. numpy is used for matrix multiplication.
2. inv is required for solving inverse of a matrix


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

# Defining the functions of coupled-system
def f1(y1_o,y1_n,y2_n,y3_n,h):
	return (y1_o + h*(-0.04*y1_n + y2_n*y3_n*1e4) - y1_n)

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

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

# Defining different parameters of Newton-Raphson method i.e. function representing Jacobian,Ynew and Yguess
def jac(y1_o,y2_o,y3_o,y1_n,y2_n,y3_n,h):


	# Jacobian matrix for 3 eqn\'s and 3 varaibles a 3x3 matrix
	J = np.ones((3,3))
	Δy = 1e-8    # Step size

	# Row 1
	J[0,0] = (f1(y1_o,y1_n+Δy,y2_n,y3_n,h)-f1(y1_o,y1_n,y2_n,y3_n,h))/Δy 
	J[0,1] = (f1(y1_o,y1_n,y2_n+Δy,y3_n,h)-f1(y1_o,y1_n,y2_n,y3_n,h))/Δy 
	J[0,2] = (f1(y1_o,y1_n,y2_n,y3_n+Δy,h)-f1(y1_o,y1_n,y2_n,y3_n,h))/Δy 

    # Row 2
	J[1,0] = (f2(y2_o,y1_n+Δy,y2_n,y3_n,h)-f2(y2_o,y1_n,y2_n,y3_n,h))/Δy 
	J[1,1] = (f2(y2_o,y1_n,y2_n+Δy,y3_n,h)-f2(y2_o,y1_n,y2_n,y3_n,h))/Δy 
	J[1,2] = (f2(y2_o,y1_n,y2_n,y3_n+Δy,h)-f2(y2_o,y1_n,y2_n,y3_n,h))/Δy

    # Row 3
	J[2,0] = (f3(y3_o,y1_n+Δy,y2_n,y3_n,h)-f3(y3_o,y1_n,y2_n,y3_n,h))/Δy 
	J[2,1] = (f3(y3_o,y1_n,y2_n+Δy,y3_n,h)-f3(y3_o,y1_n,y2_n,y3_n,h))/Δy
	J[2,2] = (f3(y3_o,y1_n,y2_n,y3_n+Δy,h)-f3(y3_o,y1_n,y2_n,y3_n,h))/Δy

	return J


# Defining the guess values to proceed
y1_o = 1
y2_o = 0
y3_o = 0

# Defining the Newton-Raphson solver parameters
err = 1e9
alpha = 1
tol = 1e-12
t_start=0
t_stop=600
n=1500
h=(t_stop-t_start)/(n-1)
t=np.arange(t_start,t_stop,h)
L_1 = []
L_2 = []
L_3 = []
repitation = []
iter =1

#Assigning value to proceed
Y_o = np.ones((3,1))
Y_o[0] = y1_o
Y_o[1] = y2_o
Y_o[2] = y3_o

#Now attributing these values for F
F = np.copy(Y_o)

# Initiating Newton raphson method loop
for i in range(0,len(t)):
	# Inner iterative solver loop
		while (err>tol):
			J = jac(y1_o,y2_o,y3_o,Y_o[0], Y_o[1], Y_o[2],h)

			F[0] = f1(y1_o,Y_o[0],Y_o[1],Y_o[2],h)
			F[1] = f2(y2_o,Y_o[0],Y_o[1],Y_o[2],h)
			F[2] = f3(y3_o,Y_o[0],Y_o[1],Y_o[2],h)
			
			Y_n = Y_o - alpha*np.matmul(inv(J), F)

			err = max(abs(Y_n-Y_o))
			print(err)
			
			# Updating the guess values for a new iteration
			Y_o = Y_n
	
		
		# Prompt message for each time-step

		log_message = \'repitation No. = {0} y1_n = {1} y2_n = {2} y3_n = {3}\' .format(iter,Y_n[0],Y_n[1],Y_n[2])
		print(log_message)

		# Storing the values for plotting
		L_1.append(Y_n[0])
		L_2.append(Y_n[1])
		L_3.append(Y_n[2])

		iter = iter+1
		repitation.append(iter)
		Y_o = Y_n
		y1_o = Y_n[0]
		y2_o = Y_n[1]
		y3_o = Y_n[2]

		err = 1e9
		
#Final result

fig1,ax1 = plt.subplots()
plt.plot(repitation,L_1,color = \'g\')
plt.plot(repitation,L_2,color = \'m\')
plt.plot(repitation,L_3,color = \'r\')
plt.grid(True)
xticks = [0,100,200,300,400,500,600,700,800,900,1000]
ax1.set_xticklabels(xticks,rotation=0,fontsize=12)
yticks = [0, 0.2, 0.4, 0.6, 0.8, 1, 1.2]
ax1.set_yticklabels(yticks,rotation=0,fontsize=12)
plt.xlabel(\'Time (sec)\')
plt.ylabel(\'Y\')
plt.title(\'Variation of Y with time-dt = \'+str(h),fontsize=14,fontweight=\'bold\')
plt.legend([\'y1\',\'y2\',\'y3\'],fontsize=16)
plt.show()

5. Plots ( At dt=10,1,0.1,0.01)

We observe that for above simulation time of 10min the solution obtained changes until the step size of 0.1 and remains constant after that hence dt=0.1 is considered as and opitimal step size.

The solution values are given below as y1,y2 and y3 as :-