Multivariate Newton Rhapson Solver

INTRODUCTION:

The Newton-Raphson method discussed above for solving a single-variable equation can be generalized for solving multivariate equation systems. A stiff ODE is the one where we use implicit methods to integrate the particular system consisting of non linear equations. The particular ODE\'s shown below replicates a reacting system where we solve for each and every species:

MATHEMATICAL FORMULATION:-

1) INTEGRATING USING IMPLICIT/ BACKWARD DIFFERENCING METHOD:

Let us integrate the first equation of the coupled non linear system making use of the implicit / backward differencing method and also rearrange the equation:

2) SETTING UP A COLUMN MATRIX:

The first equation when we solve relates to a vector and we get a set of vectors for each and every time step. The function is evaluated for the current time step. For a time step value n = 1 , we get the following vectors:

All of these elements are now packaged into a single column vector F given by:

3) MULTIVARIATE NEWTON RAPHSON TECHNIQUE:

We have multiple equations and we make use of multivariate newton raphson technique. We start with an initial array of values and we put them in a column matrix. The formula is given below:

α - Relaxation Factor

F - function (column vector)

J - Jacobian Matrix

Representing all the values as a column matrix we get the following:

• All possible derivatives for the matrix are written making use of the Jacobian
• For y initially we start with guess values and then get new values for y.
• The new values of y are then again stored as initial guess values and the process continues to the next iteration and so on.

4) NUMERICAL DIFFERENTIATION OF THE JACOBIAN MATRIX:

For differentiating our mathematical function numerically inside the Jacobian Matrix we use forward difference formula. For a function f1 the formula is provided as shown below:

Similiarly we differentiate for all possible derivatives within the Jacobian Matrix. The forward differencing is only possible for a small incremental value of dy. The value must be made as small as possible to provide better results

CODE IMPLEMENTATION IN PYTHON:

\"\"\"
Program to solve system of coupled non linear equations
		1.		dy1/dt = -0.04*y1 + (10^4)*y2*y3
		2.		dy2/dt = 0.04*y1 -(10^4)*y2*y3 -(3*10^7)**y2
		3.		dy3/dt = 3*(10^7)**y2
\"\"\"

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

#System of non linear equations


def f1( y1_old, y1_new, y2_new, y3_new, dt ):


	\"\"\"
	First Non Linear Equation
	\"\"\"
	return (y1_old + dt*(-0.04*y1_new+pow(10,4)*y2_new*y3_new) - y1_new)



def f2( y2_old, y1_new, y2_new, y3_new, dt ):


	\"\"\"
	Second Non Linear Equation
	\"\"\"
	return (y2_old+dt*(0.04*y1_new-(pow(10,4)*y2_new*y3_new)-(3*pow(10,7)*pow(y2_new,2)))-y2_new)


def f3( y3_old, y1_new, y2_new, y3_new, dt ):


	\"\"\"
	Third Non Linear Equation
	\"\"\" 
	return (y3_old+dt*(3*pow(10,7)*pow(y2_new,2))-y3_new)


#Jacobian Matrix having all possible derivatives for our Non Linear ODE system


def jacobian( y1_old, y2_old, y3_old, y1_new, y2_new, y3_new, dt ):
	

	\"\"\"
	Define the Jacobian
	\"\"\"
	J = np.ones((3,3))

	#Step size
	dy = 1e-6


	#3 by 3 unit matrix 
	#Numerical differentiation of the Jacobian Matrix
	#Row 1
	J[0, 0] = (f1( y1_old, y1_new+dy, y2_new, y3_new, dt ) - f1( y1_old, y1_new, y2_new, y3_new, dt )) / dy
	J[0, 1] = (f1( y1_old, y1_new, y2_new+dy, y3_new, dt ) - f1( y1_old, y1_new, y2_new, y3_new, dt )) / dy
	J[0, 2] = (f1( y1_old, y1_new, y2_new, y3_new+dy, dt ) - f1( y1_old, y1_new, y2_new, y3_new, dt )) / dy

	#Row 2
	J[1,0] = (f2( y2_old, y1_new+dy, y2_new, y3_new, dt ) - f2( y2_old, y1_new, y2_new, y3_new, dt )) / dy
	J[1,1] = (f2( y2_old, y1_new, y2_new+dy, y3_new, dt ) - f2( y2_old, y1_new, y2_new, y3_new, dt )) / dy
	J[1,2] = (f2( y2_old, y1_new, y2_new, y3_new+dy, dt ) - f2( y2_old, y1_new, y2_new, y3_new, dt )) / dy

	#Row 3
	J[2,0] = (f3( y3_old, y1_new+dy, y2_new, y3_new, dt ) - f3( y3_old, y1_new, y2_new, y3_new, dt )) / dy
	J[2,1] = (f3( y3_old, y1_new, y2_new+dy, y3_new, dt ) - f3( y3_old, y1_new, y2_new, y3_new, dt )) / dy 
	J[2,2] = (f3( y3_old, y1_new, y2_new, y3_new+dy, dt ) - f3( y3_old, y1_new, y2_new, y3_new, dt )) / dy

	return J

# Providing initial guess values for performing Newton iterations


y1_old = 1

y2_old = 0

y3_old = 0

Error = 9e9

#Numpy column vector

Y_old = np.ones((3,1))
Y_old[0] = 1
Y_old[1] = 0
Y_old[2] = 0


#Numpy column vector
Y_new = np.ones((3,1))
Y_new[0] = 1
Y_new[1] = 0
Y_new[2] = 0

F = np.copy(Y_old)


#Defining Newton Raphson parameters

alpha = 1
tol = 1e-6
t = 600
dt = 0.1
n = int(t/dt)

# Setting up the plot 
Species_1 = []
Species_2 = []
Species_3 = []
iteration = []

iter = 0.1



for i in range(0,n):
	
	while Error > tol:

		#Computing Jacobian old values are used

		J = jacobian( y1_old, y2_old, y3_old, Y_old[0], Y_old[1], Y_old[2], dt)


		#Computing the F vector old values are used

		F[0] = f1(y1_old,Y_old[0], Y_old[1], Y_old[2],dt)
		F[1] = f2(y2_old,Y_old[0], Y_old[1], Y_old[2],dt)
		F[2] = f3(y3_old,Y_old[0], Y_old[1], Y_old[2],dt) 


		#Computing New values

		Y_new = Y_old - alpha*(np.matmul(inv(J),F))

		#Computing Maximum Absolute Errors

		Error = np.max(np.abs(Y_new-Y_old))
		print (\'Error =\',Error)


		#Updating old values 

		Y_old = Y_new	

	#Final Results stored as guess values to perform repeated iterations
	Y_old = Y_new	
	y1_old = Y_new[0]
	y2_old = Y_new[1]
	y3_old = Y_new[2]  
	Error = 9e9

		


	#log message
	log_message = \'iteration #= {0} y1= {1} y2= {2} y3= {3}\'.format(iter, Y_new[0], Y_new[1], Y_new[2])
	print(log_message)

	Species_1.append(Y_new[0])
	Species_2.append(Y_new[1])
	Species_3.append(Y_new[2])

	iter = iter +1
	iteration.append(iter)
	
		

# Setting up the plot for our final species

plt.plot(iteration,Species_1,color = \'black\')

plt.plot(iteration,Species_2,color = \'yellow\')

plt.plot(iteration,Species_3,color = \'purple\')

plt.xlabel(\'Number of cycles TIME (secs)\')

plt.ylabel(\'Initial Guess values for each Species (Y)\')

plt.grid (\'on\')

plt.legend([\'Species 1 (Y1)\', \'Species 2 (Y2)\', \'Species 3 (Y3)\'])

plt.title(\'Multivariate Newton Raphson Method for dt = 0.1\')

plt.show()

OUTPUT:

• Making use of the implicit/ backward differencing scheme we see that our solution doesn\'t blow up for longer cycles of time and widely remains stable.
• Keeping the value of dt < 0.1 only the cycles kepping adding up but the solution remains stable and there is no change in the species values.
• Reducing the time step increases the number of cycles but there will be a drawback on the amount of time we spend on running the simulations. Most important thing is to spend less time for computation for better results

• We see that as the value of dt keeps increasing the guess values seem to not align with the initial guess values which we assigned for the Newton Raphson problem. There seems to be a deviation in the plot and when we further increase the value of dt the species will vanish completely. Also less time is spent in the computation process but it lacks the accuracy in results.

CONCLUSION:

• More iterations needs to be run for better accuracy in the results for simulation. Having the value of dt< 0.1 provides accurate results for our simulation but there is a necessity for bringing down the time required for computation.