Multivariate Newton-Raphson Solver for a Coupled Non-linear equations

Multivariate Newton Raphson Solver:

Objective: Writing a code of Multivariate Newton Raphson solver for the given coupled non-linear equations. The equations are given in the form that we cannot solve the variables y1, y2, and y3 separately and this means the word coupled. In these equations, variables are multiplied to itself or multiplied by other variables, this means the non-linear. Follow the given initial values for the variables for the further process. Take alpha as 1.

Multivariate Newton Raphson Solver

• The most widely used method for root-finding is the Newton-Raphson method.
• By starting with the initial guess x_i, a tangent can be extended from the point `[x_i,f(x_i)]`. The point where this tangent crosses the x-axis is the best estimate of the root.
• Solving the combustion process, reaction mechanisms and some other processes which contain a large number of data or equations are not that simple to solve with hand, and not easy with the code too.
• But including a mathematical model to a particular code will do the job easy for us.
• The coupled equations in the real sense, let us a combustion process is happening, the species producing from that process are not so much different from other species. But each species definitely depends on the other species.
• If we take that data as an equation, those equations are seemed to be coupled.
• Handling that data can be done by dong integration. but it\'ll take a huge amount of time. To avoid this, we are introducing the numerical methods in these types of cases.
• The Multivariate Newton Raphson Method is one of those methods which gives precisely accurate results.
• The Multivariate Newton Raphson Method is a root-finding algorithm that produces successively better approximations to the roots (or zeroes) of a real-valued function. It actually came from the Univariate Newton Raphson Method as:
• The roots are going to find by using the iterative techniques, here we use the implicit scheme.

`y_((n)ew) = y_(guess) - alpha*[(f(y_(guess)))/(f\'(y_(guess)))]`

Where y_new = new iteration value of the corresponding variable

y_guess = the value that came from the last iteration of a corresponding variable

alpha = over-relaxation factor

• Now, as the univariate method, the multivariate method is the same with some extra features included.

`y_((n)ew) = y_(guess) - alpha*(J^(-1))*F`

Where y_new = new iteration value of the corresponding variable

y_guess = the value that came from the last iteration of a corresponding variable

alpha = over-relaxation factor

J = jacobian matrix of the coupled equation

• The given coupled non-linear equations to be solved are:

`(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`

• The given initial conditions for the three variables are:

`y_1[0] = 1`

`y_2[0] = 0`

`y_3[0] = 0`

• The actual equations looks like when the differentiation applied

`(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*10^7(y_2^n)^2`

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

• The actual process of Multivariate Newton Raphson method is finding the roots of zeros.
• Now, we need to rearrange the above equations to a system of root finding functions.
• Applying the Backward differencing method for the ODE\'s.

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

`f_2^n = y_2^n - Deltat(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 - Deltat(3*10^7(y_2^n)^2) - y_3^(n-1) = 0`

Where,

`y_n = y^(New)`

`y^(n-1) = y^(old)`

• After setting up the equations for all the three functions, now we need to find the jacobian matrix using these matrices.

`J = [[(partialf_1^n)/(partialy_1^n),(partialf_1^n)/(partialy_2^n),(partialf_1^n)/(partialy_3^n)],[(partialf_2^n)/(partialy_1^n),(partialf_2^n)/(partialy_2^n),(partialf_2^n)/(partialy_3^n)],[(partialf_3^n)/(partialy_1^n),(partialf_3^n)/(partialy_2^n),(partialf_3^n)/(partialy_3^n)]]`

• The final equation or the general equation for the Multivariate Newton Raphson method is written as:

`y_(N+1)^n = y_(N)^n - alpha*[(J^(-1))*F]^n`

Where F is a functional or final matrix

Y_N+1 is a new valued matrix

Y_N is a old valued matrix

`F = [[f_1],[f_2],[f_3]]`

N+1 = New iteration

N = Old iteration

`alpha` = Relaxation factor

Code:

• The parameters used for the solver are

Defining a function as = `f_i = (partialf_i)/(partialy_i) = (f_i(y_i+Deltay_i) - f_i(y_i))/(Deltay_i)`

`Deltay_i = 10^-5`, in code `Deltay` is represented as dy.

Here i = 1,2, and 3

`Delta t = 0.5`, in code `Delta t` is represented as dt.

Simulation time = 600 seconds

Error = 9e9

tolerance = 1e9

Relaxation factor `alpha` = 1

\"\"\"
A python code to solve the given system of coupled non-linear equations using Multivariate Newton-Raphson Method
Coupled equations are
dy1/dt = -0.04*y1 + 10^4*y2*y3 --> eq.1
dy2/dt = 0.04*y1 - 10^4*y2*y3 - 3*10^7*y2^2 --> eq.2
dy3/dt = 3*10^7*y2^2 --> eq.3
\"\"\"

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

# time interval data
t_start = 0
t_end = 600
n = 1200
# t_end = 10*60 (seconds)
dt = (t_end-t_start) / n # time step

# initial guess values of co-efficients
y1 = 1
y2 = 0
y3 = 0

def f1 (y1_old, y2_old, y3_old, y1, y2, y3, dt):
	return y1_old + (-0.04*y1 + pow(10,4)*y2*y3)*dt - y1

def f2 (y1_old, y2_old, y3_old, y1, y2, y3, dt):
	return y2_old + (0.04*y1 - pow(10,4)*y2*y3 - 3*pow(10,7)*pow(y2,2))*dt - y2

def f3 (y1_old, y2_old, y3_old, y1, y2, y3, dt):
	return y3_old + (3*pow(10,7)*pow(y2,2))*dt - y3

def jac (y1_old, y2_old, y3_old, y1, y2, y3, dt):
	j = np.zeros((3,3))
	dy = 1e-5

	# row - 1
	j[0,0] = (f1(y1_old, y2_old, y3_old, y1+dy, y2, y3, dt) - f1(y1_old, y2_old, y3_old, y1, y2, y3, dt))/dy
	j[0,1] = (f1(y1_old, y2_old, y3_old, y1, y2+dy, y3, dt) - f1(y1_old, y2_old, y3_old, y1, y2, y3, dt))/dy
	j[0,2] = (f1(y1_old, y2_old, y3_old, y1, y2, y3+dy, dt) - f1(y1_old, y2_old, y3_old, y1, y2, y3, dt))/dy

	# row - 2
	j[1,0] = (f2(y1_old, y2_old, y3_old, y1+dy, y2, y3, dt) - f2(y1_old, y2_old, y3_old, y1, y2, y3, dt))/dy
	j[1,1] = (f2(y1_old, y2_old, y3_old, y1, y2+dy, y3, dt) - f2(y1_old, y2_old, y3_old, y1, y2, y3, dt))/dy
	j[1,2] = (f2(y1_old, y2_old, y3_old, y1, y2, y3+dy, dt) - f2(y1_old, y2_old, y3_old, y1, y2, y3, dt))/dy

	# row - 3
	j[2,0] = (f3(y1_old, y2_old, y3_old, y1+dy, y2, y3, dt) - f3(y1_old, y2_old, y3_old, y1, y2, y3, dt))/dy
	j[2,1] = (f3(y1_old, y2_old, y3_old, y1, y2+dy, y3, dt) - f3(y1_old, y2_old, y3_old, y1, y2, y3, dt))/dy
	j[2,2] = (f3(y1_old, y2_old, y3_old, y1, y2, y3+dy, dt) - f3(y1_old, y2_old, y3_old, y1, y2, y3, dt))/dy

	return j

def Newton_Raphson(y1, y2, y3, y1_guess, y2_guess, y3_guess, dt):
	# calculating the numerical jacobian in a vector form and updating the old y values by indicating as \'Y\'

	Y_old = np.zeros((3,1))
	Y_old[0] = y1_guess 
	Y_old[1] = y2_guess
	Y_old[2] = y3_guess

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

	# definig the solver parameters
	error = 9e9
	tol = 1e9
	alpha = 1

	iter = 1

	while error >= tol:
		# defining how jacobi sets for the equation
		j = jac(y1, y2, y3, Y_old[0], Y_old[1], Y_old[2], dt)

		F[0] = f1(y1, y2, y3, Y_old[0], Y_old[1], Y_old[2], dt)
		F[1] = f2(y1, y2, y3, Y_old[0], Y_old[1], Y_old[2], dt)
		F[2] = f3(y1, y2, y3, Y_old[0], Y_old[1], Y_old[2], dt)

		# After guessing the y values, now computing the new y values
		Y_new = Y_old - alpha*(np.matmul(inv(j), F))

		# taking the maxing possible error from the vector matrix of Y
		error = np.max(np.abs(Y_new-Y_old))

		# equating old values as new values(updating)
		Y_old = Y_new

		# log_message for the printing of number of iterations
		log_message = \'iteration = {0} y1 = {1} y2 = {2} y3 = {3}\'.format(iter, Y_new[0], Y_new[1], Y_new[2])
		#log_message = \'Time = {0} sec, y1 = {1}, y2 = {2}, y3 = {3}\'.format(round(t[k], 1),round(Y_new[0][0],3),round(Y_new[1][0],3),round(Y_new[2][0],3)
		print (log_message)
		iter = iter + 1

	return [Y_new[0], Y_new[1], Y_new[2]]

def implicit_euler (y1, y2, y3, t_start, t_end, dt):

	t = np.arange(t_start, t_end, dt)
	time = len(t)

	Y_1 = np.zeros(time)
	Y_2 = np.zeros(time)
	Y_3 = np.zeros(time)

	Y_1[0] = y1
	Y_2[0] = y2
	Y_3[0] = y3

	y1_guess = 0
	y2_guess = 0
	y3_guess = 0

	for i in range(1, time):
		Y_1[i], Y_2[i], Y_3[i] = Newton_Raphson(Y_1[i-1], Y_2[i-1], Y_3[i-1], y1_guess, y2_guess, y3_guess, dt)

		# final output result

		y1_guess = Y_1[i]
		y2_guess = Y_2[i]
		y3_guess = Y_3[i]

	return [t, Y_1, Y_2, Y_3]

t, y1, y2, y3 = implicit_euler(1,0,0,0,600,0.05)

plt.plot (t, y1, color = \'blue\')
plt.plot (t, y2, color = \'green\')
plt.plot (t, y3, color = \'magenta\')
plt.grid (\'on\')
plt.xlabel (\'Simulation Time (sec)\')
plt.ylabel (\'y_i\')
plt.legend([\'y1\',\'y2\',\'y3\'])
plt.title(\'Multivariate Newton-Raphson Method for dt = 0.5\')
plt.show() 

Plot:

for dt = 0.5 (n = 1200)

for dt = 0.05 (n = 12000)

for dt = 5 (n = 120)

Conclusion:

• From the above graphs, we can conclude that the approximate values of y1, y2, and y3 are changing rapidly with respect to the time in the X-axis.
• These values are changing rapidly as the solution marches forward in time. This nature is observed because of the stiff character of the system of equation.

Choosing a suitable time step size:

• Firstly the time step dt is set to 0.5 then run the code and obtain a solution. For the second time changed the time step value dt to 0.05 and run the simulation. But there is no difference in the solution and plot. For the third time, the time step increased aggressively to dt = 5 and run the simulation. But there is no change at this attempt also. The simulation run times also didn\'t change much as well.
• So, finally, it says that the change in time step remains the solution same, stable and converging finally giving meaning full results.
• So, there is no effect of the time step size on the simulation.
• The reason for the simulation unconditionally stable is, we are using the implicit scheme as an iterative solver and it has no effect of time step size on it.
• Actually, the change in time step size may affect the solution if we use an explicit scheme as the iterative solver.
• Even though the implicit schemes are somewhat difficult to code, but can work extremely good st complex situations like the stiff characteristic of the non-linear coupled equations, situations like a sudden rise in temperature, pressure, and mass flow rate, etc.