Multivariate Newton Rhapson Solver

INTRODUCTION & BACKGROUND:

The given equations represent a reacting flow mechanism. Usually the ODE\'s of a reacting flow are stiff in nature, meaning the solution of the ODE varies rapidly and hence there is a need to choose `trianglet` in such a way that the solution does not go unstable. This is the main problem that is encountered when solving the ODE explicitly.

Hence, in this project the given system of ODE is to be solved using Implicit Backward Differencing Scheme with Newton Raphson convergence. The main advantage of Implicit scheme is that the solution does not become unstable at larger values of `triangle t` but there is a loss in the accuracy of the solution.

The Implicit BDF scheme is given by:

`y^(n)underset(i)() = y^(n-1)underset(i)() + h*f^n`

where \'f\' is the function, \'h\' is the timestep and \'n\' represents the current time iteration. Hence \'n-1\' term becomes the old value and other terms are in the fucntion are calculated at current time.

The Newton Raphson convergence formula is given by:

`y^(n ew) = y^(old) + alphaJ\'F`

where y represents the function at new and old timesteps, `alpha` is the under-relaxation factor, J\' is the transpose of the Jacobian of the functions and F is the function computed at the guess values.

In this project, the jacobian is calculated numerically using the slope formula: `frac{f(x+triangleh) - f(x)}{triangleh}`. Also, the effect of larger timestep is shown as well.

PYTHON CODE:

\"\"\"

THE EQUATIONS ARE:

3*x1 - cos(x2*x3) - 3/2 = 0
4*x1^2 - 625*x2^2 + 2*x3 - 1 = 0
20*x3 + e^(-x1*x2) + 9 = 0

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

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

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

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



def jacobian(y1_old, y2_old, y3_old, y1, y2, y3, dt):

	J = np.ones((3,3))

	J[0,0] = (f1(y1_old, y1+1e-6, y2, y3, dt) - f1(y1_old, y1, y2, y3, dt))/1e-6
	
	J[0,1] = (f1(y1_old, y1, y2+1e-6, y3, dt) - f1(y1_old, y1, y2, y3, dt))/1e-6

	J[0,2] = (f1(y1_old, y1, y2, y3+1e-6, dt) - f1(y1_old, y1, y2, y3, dt))/1e-6

	J[1,0] = (f2(y1_old, y1+1e-6, y2, y3, dt) - f2(y1_old, y1, y2, y3, dt))/1e-6

	J[1,1] = (f2(y2_old, y1, y2+1e-6, y3, dt) - f2(y2_old, y1, y2, y3, dt))/1e-6

	J[1,2] = (f2(y2_old, y1, y2, y3+1e-6, dt) - f2(y2_old, y1, y2, y3, dt))/1e-6

	J[2,0] = (f3(y3_old, y1+1e-6, y2, y3, dt) - f3(y3_old, y1, y2, y3, dt))/1e-6

	J[2,1] = (f3(y3_old, y1, y2+1e-6, y3, dt) - f3(y3_old, y1, y2, y3, dt))/1e-6

	J[2,2] = (f3(y3_old, y1, y2, y3+1e-6, dt) - f3(y3_old, y1, y2, y3, dt))/1e-6


	return J





t_start = 0;
t_end = 60*10;
n = 100		#No of time steps
dt = (t_end - t_start)/n-1
iter = 1

y_guess = np.ones((3,1))
y_old = np.ones((3,1))
y_new = np.ones((3,1))


#numpy col vector for F
F = np.copy(y_guess)

alfa = 1
plot_ctr = np.linspace(0,100,n)
plot_y1 = [None]*100
plot_y2 = [None]*100
plot_y3 = [None]*100
tol = 1e-8

#Initialize loop value
y1 = 1;
y2 = 0;
y3 = 0;

while iter <= n:								#time loop

	iter_raph = 1


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

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

	error = 9e9									#Initializing/resetting error before start of every conv. loop
	
	while error > tol:							#NR Convergence loop
		
		J = jacobian(y_old[0], y_old[1], y_old[2], y_guess[0], y_guess[1], y_guess[2], dt)
		#Computing F vector 
		F[0] = f1(y_old[0], y_guess[0], y_guess[1], y_guess[2], dt)
		F[1] = f2(y_old[1], y_guess[0], y_guess[1], y_guess[2], dt)
		F[2] = f3(y_old[2], y_guess[0], y_guess[1], y_guess[2], dt)

		y_new = y_guess - alfa*np.matmul(inv(J),F)

		error = np.max(np.abs(y_new - y_guess))

		y_guess = y_new

		#log_msg
		log_message = \'iteration # = {0} y1 = {1} y2 = {2} y3 = {3} error = {4}\'.format(iter_raph, y_guess[0], y_guess[1], y_guess[2], error)
		print(log_message)
		iter_raph +=1

	
	

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


	
	plot_y1[iter-1] = y_new[0]
	plot_y2[iter-1] = y_new[1]
	plot_y3[iter-1] = y_new[2]
	iter +=1
	


plt.plot(plot_ctr, plot_y1, \'-\', color = \'red\', label = \'Y1\') 
plt.plot(plot_ctr, plot_y2, \'-\', color = \'blue\', label = \'Y2\')
plt.plot(plot_ctr, plot_y3, \'-\', color = \'green\', label = \'Y3\')
plt.xlabel(\'Time (seconds)\')
plt.ylabel(\'Values\')
plt.legend()
plt.show()

SOLUTION OBTAINED:

Solution at dt = 1.2 Solution at dt = 6

Solution at dt = 60

INFERENCE:

We can clearly see that the solution is becoming less acccurate as the timestep increases in magnitude. However, the solution does not show signs of oscillations or blow up even with higher timesteps. This is advantage of Implicit BDF scheme.

For this project, a timestep of 1.2 is good enough to give fairly accurate results; On increasing the size of steps, the accuracy decreases, and on decreasing the size of steps, the simulation time and strain on computer increases.