Week 6 - Multivariate Newton Rhapson Solver

Aim

To solve functions in Multivariate Newton Rhapson sover with backward differencing method.

Multivariate Newton Rhapson

The Multivariate Newton Rhapson is sloved by using the following formula.

`Y_(n ew) = Y_(old) - alpha J^-1 F`

where

Y_new is the incremented value of the each iteration of the time loop is stored.

Y_old is the initial values are stored .

`alpha`is the relaxtion factor

`J^-1`is inverse of the Jacobian matrix

F is the function matrix.

""""
The functions for the program are
dy1/dt = -0.04y1 + 10^4*y2*y3
dy2/dt = 0.04y1 - 10^4*y2*y3 - 3*10^7*y2*y2
dy3/dt = 3*10^7*y2*y2
"""

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

#Defining the Functions in Implicit method,where h is the Time step
def f1(y1_old, y1, y2, y3, h):
	return y1_old + h*(-0.04*y1 + 1e4*y2*y3) - y1

def f2(y2_old,y1,y2,y3,h):
	return y2_old + h*(0.04*y1 - 1e4 *y2*y3 - 3e7*y2*y2) - y2
def f3(y3_old,y1,y2,y3,h):
	return y3_old + h*(3e7*y2*y2) - y3		

#Defining the Jacobian Matrix for each for the functions
def jacobian(y1_old,y2_old,y3_old,y1,y2,y3,h):

	J= np.ones((3,3))
	d =1e-6
	
	#Row 1
	J[0,0] = (f1(y1_old,y1+d,y2,y3,h) - f1(y1_old,y1,y2,y3,h))/d
	J[0,1] = (f1(y1_old,y1,y2+d,y3,h) - f1(y1_old,y1,y2,y3,h))/d
	J[0,2] = (f1(y1_old,y1,y2,y3+d,h) - f1(y1_old,y1,y2,y3,h))/d

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

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

	return J

#Initial Values
y1= 1
y2= 0
y3= 0	

Y_old = np.ones((3,1))
Y_guess = np.ones((3,1))

#Time Steps
h = 1
t_initial = 0 
t_final = 600 + h
t = np.arange(t_initial,t_final,h)
total_iter = len(t)


F= np.copy(Y_guess)
time_iter = 0

Y = np.ones((3,total_iter))
Y[0,0] = y1
Y[1,0] = y2
Y[2,0] = y3


log_message = 'flow time={0}s,iteration ={1},y1={2},y2={3},y3={4}'.format(time_iter,'initial',y1,y2,y3)
print(log_message)

#Outer Time Loop
for i in range(1,total_iter):
	Y_old[0] = y1
	Y_old[1] = y2
	Y_old[2] = y3

	Y_guess[0] = y1
	Y_guess[1] = y2
	Y_guess[2] = y3

	error  =9e9
	alpha = 1
	tol = 1e-8

	iter =1
	while error>tol:
		#computing the Jacobian 
		J = jacobian(Y_old[0],Y_old[1],Y_old[2],Y_guess[0],Y_guess[1],Y_guess[2],h)

		F[0] = f1(Y_old[0],Y_guess[0],Y_guess[1],Y_guess[2],h)
		F[1] = f2(Y_old[1],Y_guess[0],Y_guess[1],Y_guess[2],h)
		F[2] = f3(Y_old[2],Y_guess[0],Y_guess[1],Y_guess[2],h)

		#Computing the Values
		Y_new = Y_guess - alpha*np.matmul(inv(J),F)

		#Error 
		error = np.max(np.abs(Y_guess-Y_new))

		#Updating the Values
		Y_guess = Y_new
		iter = iter+1

	y1=Y_new[0]
	y2=Y_new[1]
	y3=Y_new[2]	

	time_iter = time_iter+1
	Y[0,time_iter] = Y_new[0]
	Y[1,time_iter] = Y_new[1]
	Y[2,time_iter] = Y_new[2]


	log_message = 'flow time={0}s,iteration ={1},y1={2},y2={3},y3={4}'.format(time_iter*h,iter-1,Y_new[0],Y_new[1],Y_new[2])
	print(log_message)

#plotting the funcions with respect to time 
plt.plot(t,Y[0,:],color ='red')
plt.plot(t,Y[1,:],color ='blue')
plt.plot(t,Y[2,:],color ='orange')
plt.title('Multivariate Newton Rhapson')
plt.xlabel('Flow Time(seconds)')
plt.ylabel('Functions')
plt.legend(['Y1','Y2','Y3'])
plt.show()

Time Step	y1 at 600s	y2 at 600s	y3 at 600s
0.01	0.39992672	2.63125125e-06	0.60007065
0.1	0.39997578	2.63177957e-06	0.60002159
1	0.40045615	2.63695669e-06	0.59954121
10	0.40476791	2.68377845e-06	0.59522941

When the time step of 1 is changed to 0.01,0.1,10 the computational time has increased with respect to the accuracy.If the purpose of the computation is with the accutacy then decreasing the timestep is the valid but if the compuatuon is for the faster time then increasing the timestep is valid with loss in accuracy. For this the Time step 1 is ideal giving ideal accuracy with ideal computational time.