Week 6 - Multivariate Newton Rhapson Solver

Multivariate Newton Rapshon Method:-

• In numerical analysis, Newton's Method also known as the Newton-Raphson method is a root.
• Finding algorithms which produce successively better approximation to the root or zeros of a real values function.
• The most basic version start with a single variable function defined for a real variable x, the function f defined for a real variable x, the function f' and an initial guess xoxo for the root of f
• The idea is to start with an initial guess which is reasonably close to the true root, then to approximate the function by it tangent line using calculus and finally to compute the x-intercept of its tangent line by elementary algebra.

• This x-intercept will typically be a better approximation to the original function root than the first guess, and the method can be iterated.

Univariate Newton Raphson Method

y_new = y_guess - alpha*[f(y_guess)/f'(y_guess)]

y_new = yguess−α⋅(J−1⋅F)yguess-α⋅(J-1⋅F)

where y_new = new iteration value that from a corresponding variable

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

alpha = relaxation factor

J = Jacobian Matrix

Given coupled non-linear equation

dy1dt=−0.04y1+104y2y3dy1dt=-0.04y1+104y2y3

dy2dt=0.04y1+104y2y3−3⋅107⋅y22dy2dt=0.04y1+104y2y3-3⋅107⋅y22

dy3dt=3⋅107⋅y22dy3dt=3⋅107⋅y22

Initial condition

y1[0]=1,y2[0]=0,y3[0]=0y1[0]=1,y2[0]=0,y3[0]=0

Applying backward differentiation

yn1−yn−11△t=−0.04yn1+104⋅yn2.yn3y1n-y1n-1△t=-0.04y1n+104⋅y2n.y3n

yn2−yn−12=0.04yn1−104yn2yn3−3⋅107(yn2)2y2n-y2n-1=0.04y1n-104y2ny3n-3⋅107(y2n)2

dy3dt=3⋅107⋅(yn2)2dy3dt=3⋅107⋅(y2n)2

Mutlivarient Newton Rapshon Method

fn1=yn1+△t(0.04yn1−104yn2yn3)−yn−11=0f1n=y1n+△t(0.04y1n-104y2ny3n)-y1n-1=0

fn2=yn2−△t(0.04yn1−104yn2yn3−3⋅107(yn2)2)−yn−12=0f2n=y2n-△t(0.04y1n-104y2ny3n-3⋅107(y2n)2)-y2n-1=0

fn3=yn3−△t(3⋅107(yn2)2)−yn−13=0f3n=y3n-△t(3⋅107(y2n)2)-y3n-1=0

To find Jacobian Matrix

Since we have three equation and three variable. Therefore, there will be three rows and columns respectively.

J=⎡⎢ ⎢ ⎢ ⎢ ⎢ ⎢⎣∂fn1∂yn1∂fn1∂yn2∂fn1∂yn3∂fn2∂yn1∂fn2∂yn2∂fn2∂yn3∂fn3∂yn1∂fn3∂yn2∂fn3∂yn3⎤⎥ ⎥ ⎥ ⎥ ⎥ ⎥⎦J=[∂f1n∂y1n∂f1n∂y2n∂f1n∂y3n∂f2n∂y1n∂f2n∂y2n∂f2n∂y3n∂f3n∂y1n∂f3n∂y2n∂f3n∂y3n]

y(new)n=ynold+α[(J−1)F]ny(new)n=yoldn+α[(J-1)F]n

Here, F=⎡⎢⎣f1f2f3⎤⎥⎦F=[f1f2f3]

Code:-

""
Program to solve the coupled non linear equation

dy1/dt = -0.04*y1 + pow(10,4)*y2*y3
dy2/dt = 0.04*y1 - pow(10,4)*y2*y3 - 3*pow(10,7)*y2*y2
dy3/dt = 3*pow(10,7)*y2*y2
"""
import numpy as np 
from numpy.linalg import inv
import matplotlib.pyplot as plt 

# Time Interval
t_start = 0
t_end = 600
n = 6000
dt = (t_end - t_start)/n

iteration = 0

# Initial Guess Value
y1 = 1
y2 = 0
y3 = 0

# Equation
def f1(y1_old,y2_old,y3_old,y1,y2,y3,dt):
	""" First Non - Linear Equation"""
	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):
	""" Second Non - Linear Equation"""
	return y2_old + (0.04*y1 - pow(10,4)*y2*y3 - 3*pow(10,7)*y2*y2)*dt - y2

def f3(y1_old,y2_old,y3_old,y1,y2,y3,dt):
	""" Third Non - Linear Equation"""
	return y3_old + (3*pow(10,7)*y2*y2)*dt - y3

# Jacobian Matrix
def jacobian(y1_old,y2_old,y3_old,y1,y2,y3,dt):
	J = np.zeros((3,3))
	h = 1e-5

	# First Row
	J[0,0] = (f1(y1_old,y2_old,y3_old,y1+h,y2,y3,dt)-f1(y1_old,y2_old,y3_old,y1,y2,y3,dt))/h
	J[0,1] = (f1(y1_old,y2_old,y3_old,y1,y2+h,y3,dt)-f1(y1_old,y2_old,y3_old,y1,y2,y3,dt))/h
	J[0,2] = (f1(y1_old,y2_old,y3_old,y1,y2,y3+h,dt)-f1(y1_old,y2_old,y3_old,y1,y2,y3,dt))/h

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

	# Third Row
	J[2,0] = (f3(y1_old,y2_old,y3_old,y1+h,y2,y3,dt)-f3(y1_old,y2_old,y3_old,y1,y2,y3,dt))/h
	J[2,1] = (f3(y1_old,y2_old,y3_old,y1,y2+h,y3,dt)-f3(y1_old,y2_old,y3_old,y1,y2,y3,dt))/h
	J[2,2] = (f3(y1_old,y2_old,y3_old,y1,y2,y3+h,dt)-f3(y1_old,y2_old,y3_old,y1,y2,y3,dt))/h

	return J



# Newton Rhaphson Method
def Newton_Rhaphson(y1,y2,y3,y1_guess,y2_guess,y3_guess,dt):
	# Vector Equation
	global iteration
	Y_old = np.zeros((3,1))
	Y_old[0] = y1_guess
	Y_old[1] = y2_guess
	Y_old[2] = y3_guess
	#print(Y_old)

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

	# Solver Parameter
	error = 9e9
	tol = 1e-4
	alpha = 1
	#iter = 0

	while(error>=tol):
		# Jacobi set for the equation
		iteration = iteration + 1
		J = jacobian(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)

		Y_new = Y_old - alpha*(np.matmul(inv(J),F))
		error = np.max(np.abs(Y_new-Y_old))
		print(error,tol)
		Y_old = Y_new
		#print(iter)
		log_message = 'iteration = {0} y1 = {1} y2 = {2} y3 = {3}'.format(iteration,Y_new[0],Y_new[1],Y_new[2])
		print(log_message)
		#iter = iter + 1

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

# Implicit Euler

def implicit(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_Rhaphson(Y_1[i-1],Y_2[i-1],Y_3[i-1],y1_guess,y2_guess,y3_guess,dt)
		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(1,0,0,0,600,0.1)

plt.plot(t,y1)
plt.plot(t,y2)
plt.plot(t,y3)
plt.xlabel('Simulation Time (sec)')
plt.ylabel('Variable')
plt.legend(['y1','y2','y3'])
plt.title('Multivarient Newton Ranpshon Method for dt = 0.1')
plt.grid('on')
plt.show()

Results:-

Choosing a suitable time step size:-

• The simulation has been run for three different time step that 0.1,1 and 0.01 repectively.
• But there is not much change in the graph
• There is no effect of the time step on the simulation because we have implicit method and if we have used the explicit it will get effect due to time step.
• However, the implicit method is difficult to code as certain extend as compared to explicit.
• But for choosing a suitable time step size we can conclude that time step of 0.01 is suitable because the Number of iteration is more, more iteration there will be more accurancy in the solution.
• Hence, time step size of 0.01 is ideal time step size.