Week 6 - Multivariate Newton Rhapson Solver

AIM - To incorporate multivariate newton raphson solver for a time marching problem

OBJECTIVE - To create the python code that solves for transient multivariate newton raphson method and see if the implicit euler method works for varirying time step values.

THEORY -  In transient multivariate newton raphson solver we will have to create a code that will use implicit euler method or what we call as backward differencing scheme where the values of current timestep are calculated using the values of current and past time step.

For this problem we have the following set of equations

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

`dy_2/dt = ( 0.04y_1 - 10^4y_2 y_3 - 3*10^7y_2^2`

`dy_3/dt = 3*10^7 y_2^2`

to solve this set of nonlinear equation using implicit euler method we will have to convert this equation into the form that represents the values at both the current and past timesteps.

this can be done by using the annotations like `y_i^n` denotes the values for at the current time step and `y_i^(n-1)` denoted the past time step values.

from the numerical definition of the derivative using backward differencing scheme

`dy_i/dt = (y_i^n - y_i^(n-1))/(Delta t)`

so our set of equation will become

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

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

now we will be arranging the equations by taking the terms on the LHS to the RHS

`y_1^(n-1) + Delta t *(-0.04 y_1^n + 10^4 y_2^ny_3^n) - y_1^n = 0`

`y_2^(n-1) + Delta t * ( 0.04 y_1^n - 10^4 y_2^n y_3^n - 3.10^7 (y_2^n)^2) - y_2^n = 0`

`y_3^(n-1) + Delta t * ( 3* 10^7 (y_2^n)^2) - y_3^n = 0`

now we will use this equation in our code to solve using the multivariate newton raphson method

also here we have to incorporate the numerical derivative technique to calculate the jacobian in our calculations for that case we will assume that we have 2 node that are very close to each other at a distance h , keeping h very small about 1e-6 we will find the numerical derivative for the functions.

the numerical derivaitve for the first jacobian term will be given with respect to the `y_1` as shown below

`f_i' (y_1, y_2, y_3 , sf"other variables") = (f_i(y_1 + h, y_2, y_3 , sf"other variables") - f_i (y_1, y_2, y_3 , sf"other variables") )/h`

similarly we can find the other derivatives terms in the jacobian using this logic 

Now to write the code first we need to create funtions for each equation that we have also we will create a function to solver for the jacobian .Then we will be having our outer time loop inside which we will have our newton raphson loop(while loop)that will solve iterativesly untill a desired convergence is reached. finally we will be plotting the value for each variables at the end of each time step.

CODE - 

Following is the code that is used to solve as set of nonlinear coupled equations using a multivariate newton raphson method.

"""
For this project the followin set of equations will be 
solved using multivariate newton raphson method.

dy1/dt = -0.04y1 +10^4y2y3
dy2/dt = 0.04y1 - 10^4y2y3 - 3.10^7y2^2
dy3/dt = 3.10^7 y2^2

modules to be imported
1. numpy module for matrix calculations
2. inv module from the numpy.linalg for inv of matrix
"""

import numpy as np
import numpy.linalg as npla
import matplotlib.pyplot as plt

def f1(y1,y2,y3,y1_old,dt):
	"""
	first non-linear equation
	"""
	return y1_old + dt*(-0.04*y1 + pow(10,4)*y2*y3) - y1

def f2(y1,y2,y3,y2_old,dt):
	"""
	second non-linear equation
	"""
	return y2_old + dt*(0.04*y1 - pow(10,4)*y2*y3 - 3*pow(10,7)*pow(y2,2)) - y2

def f3(y1,y2,y3,y3_old,dt):
	"""
	third non-linear equation
	"""
	return y3_old + dt*(3*pow(10,7)*pow(y2,2)) - y3

def jacobian(y1,y2,y3,y1_old,y2_old,y3_old,dt,h):
	# creating the jacobian matrix
	J = np.ones((3,3))

	# to calculate jacobian numerically we will consider two nodes at distance h apart 
	
	# row 1
	J[0,0] = (f1(y1+h,y2,y3,y1_old,dt) - f1(y1,y2,y3,y1_old,dt))/h 
	J[0,1] = (f1(y1,y2+h,y3,y1_old,dt) - f1(y1,y2,y3,y1_old,dt))/h
	J[0,2] = (f1(y1,y2,y3+h,y1_old,dt) - f1(y1,y2,y3,y1_old,dt))/h

	# row 2
	J[1,0] = (f2(y1+h,y2,y3,y2_old,dt) - f2(y1,y2,y3,y2_old,dt))/h 
	J[1,1] = (f2(y1,y2+h,y3,y2_old,dt) - f2(y1,y2,y3,y2_old,dt))/h
	J[1,2] = (f2(y1,y2,y3+h,y2_old,dt) - f2(y1,y2,y3,y2_old,dt))/h

	# row 3
	J[2,0] = (f3(y1+h,y2,y3,y3_old,dt) - f3(y1,y2,y3,y3_old,dt))/h 
	J[2,1] = (f3(y1,y2+h,y3,y3_old,dt) - f3(y1,y2,y3,y3_old,dt))/h
	J[2,2] = (f3(y1,y2,y3+h,y3_old,dt) - f3(y1,y2,y3,y3_old,dt))/h

	return J


# intial guess values
y1=1
y2=0
y3=0
y1_old = 1
y2_old = 0
y3_old = 0

# numpy inital guess array Y_old
Y_old = np.ones((3,1))
Y_old[0] = y1_old
Y_old[1] = y2_old
Y_old[2] = y3_old

#numpy inital array Y
Y = np.ones((3,1))
Y[0] = y1
Y[1] = y2
Y[2] = y3

# nodal distance h
h= 1e-6

# relaxation factor
alpha = 1
# tolerance 
tol = 1e-9
# total time in seconds
total_t = 600


dt = 1e-1 # timestep value
n_t = total_t/dt # total number of timesteps

# initial array F 
F = np.copy(Y_old)

# array to store the value of the variables
Y1_final = []
Y2_final = []
Y3_final = [] 

# time-step counter
t_s = 1

#time loop
for i in range (1,(int(n_t)+1)):
	
	# iteration counter
	iter = 1
	# initial error 
	error = 9e9
	
	# newton raphson loop
	while error > tol:

		J = jacobian(Y[0],Y[1],Y[2],Y_old[0],Y_old[1],Y_old[2],dt,h)

		# computing the F vector 
		F[0] = f1(Y[0],Y[1],Y[2],Y_old[0],dt)
		F[1] = f2(Y[0],Y[1],Y[2],Y_old[1],dt)
		F[2] = f3(Y[0],Y[1],Y[2],Y_old[2],dt)

		#computing new values
		Y_new = Y - alpha*np.matmul(npla.inv(J),F)
		# computing the error
		error = np.max(np.abs(Y-Y_new))
		#updating the values of Y
		Y = Y_new
		
		# log message for each iteration
		log_message1 = "iteration # = {0} ,y1 = {1} ,y2 = {2} ,y3 = {3}".format(iter,Y_new[0],Y_new[1],Y_new[2])
		print(log_message1)
		# update iteration number
		iter = iter +1
	
	# updating the Y_old values 
	Y_old = Y_new
	# storing  the value of the new values in the array at the end of each time step
	Y1_final.append(Y_new[0]) 
	Y2_final.append(Y_new[1])
	Y3_final.append(Y_new[2])
	
	
	# log message for each time-step
	log_message = "timestep # = {0} ,y1 = {1} ,y2 = {2} ,y3 = {3}".format(t_s,Y[0],Y[1],Y[2])
	print(log_message)

	# update the timestep number
	t_s = t_s + 1

# plotting the trend of the variables
plt.plot(Y1_final)
plt.plot(Y2_final)
plt.plot(Y3_final)
plt.xlabel('number of time-steps')
plt.ylabel('variables')
plt.legend(["Y1","Y2","Y3"])
plt.title("For dt = "+str(dt)+" sec  with total time = "+str(total_t)+" seconds ")
plt.grid("on")
plt.show()

the output that we get form running this code is shown below

in the code we can vary the value of time-step 'dt' and still we will get the same result.

In our case we have used 'dt = 0.1' as this value is good enough to get good results and is also quick , we can further increase the value of dt and get result more quicker but the values will be a little approximate and not that accurate though the solution will be always stable.

Further decreasing the timestep will resulting more computation cost as well as time.

Below are some output results with varying the value of dt in the code at line number 91.

For dt = 1 

For dt = 0.5

For dt = 0.01

for dt = 0.01 it takes a lot of time to solve and not recomended as it is costlier. 

CONCLUSION -

1. Here in this project we have successfully applied multivariate newton raphson technique to solve a set of nonliner equations.
2. We have learned how to incrporate backward differencing scheme into our code as well. 
3. Using the implicit method we can solve the problem with changing the timestep value and can still get a stable result.