Multivariate Newton-Raphson method

The purpose of this assignment is to create a Python program including a Multivariate Newton Rhapson Solver, to solve a non-linear coupled differential system. This method is really useful for stiff systems, where the explicit solver are unstable.

The system of equations is:

`dy_1/dt = -0.04*y_1 + 10^4*y_2*y_3`

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

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

 

This multivariate method is based on the original Newton Rhapson solver. In each time step, the iteration process for the multivariate Newton Rhapson method takes place. For that, the following equation is used iteratively, starting  by guessing a value for Xold.

`X^(n+1) = X^n - alpha*J^-1F`

In here, J is the Jacobian matrix, and X and F are column vectors of the dimension of the system; and alpha is the relaxation factor. For this program, the Jacobian has been computed numerically, as not always the analytical derivatives are available. As an example, we leave here how to compute a Jacobian component numerically.

`(df_i)/(dx_j) = (f_i(x_j + epsilon) - f_i(x_j))/epsilon`

 

The solution for the previous system is represented in the following image.

 

"""
Program to solve a system of differential 
non-linear coupled equations

	dy1/dt = -0.04*y1 + 10^4*y2*y3
	dy2/dt = 0.04*y1 - 10^4*y2*y3 -3*10^7*y2^2
	dy3/dt = 3*10^7*y2^2

# We will use a BDF using Multivariate Newton Raphson Method

	f1 = y1(n-1) + dt*(-0.04*y1(n) + 10^4*y2(n)*y3(n)) - y1(n) = 0
	f2 = y2(n-1) + dt*(0.04*y1(n)-10^4*y2(n)*y3(n)-3*10^7*y2(n)^2) - y2(n) = 0
	f3 = y3(n-1) + dt*(3*10^7*y2(n)^2) - y3(n) = 0

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

#Problem Parameters
a1 = -0.04
a2 = 1e4
a3 = 3e7

#Variables for the time loop
dt=0.05
t_final = 10*60
t_initial = 0


def f1(y1_old,y1,y2,y3,dt):
	return y1_old + dt*(a1*y1 + a2*y2*y3) - y1

def f2(y2_old,y1,y2,y3,dt):
	return y2_old + dt*(-a1*y1-a2*y2*y3-a3*y2**2) - y2

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

def jacobian(y1,y2,y3,dt):
	J = np.ones((3,3))

	#Row 1
	J[0,0] = -0.04*dt - 1
	J[0,1] = 10**4*y3
	J[0,2] = 10**4*y2

	#Row 2
	J[1,0] = 0.04*dt 
	J[1,1] = dt*(-10**4*y3 - 6*10**7*y2) - 1
	J[1,2] = dt*(-10**4*y2)

	#Row 3
	J[2,0] = 0
	J[2,1] = dt*(6*10**7*y2)
	J[2,2] = -1

	return J

#Defining the initial and guess values
y1 = 1
y2 = 0
y3 = 0
#numpy column vector
Y_old = np.ones((3,1))
#numpy column vector
Y_guess = np.ones((3,1))

#numpy time vector
t = np.arange(t_initial,t_final,dt)
total_iter = len(t)

#numpy solution matrix
Y = np.ones((3,total_iter))

#numpy column vecotr for F
F = np.copy(Y_guess)

#Starting the Time integration Loop
time_iter = 0

for i in range(1,total_iter+1):

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

	Y_guess[0] = y1
	Y_guess[1] = y2
	Y_guess[2] = y3
	
	#Parameters for Newton Raphson solver
	error = 9e9  #Initializing error
	alpha = 1
	tol = 1e-8

	#Starting the Newton Raphson loop
	iter = 1

	while error > tol:

		#Computing the F vector - old values are used
		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)

		#Computed the Analytical Jacobian  - old values are used
		J = jacobian(Y_guess[0], Y_guess[1], Y_guess[2], dt)
		#print(J)

		# Computing the Numerical Jacobian
		n=len(Y_guess)
		eps = 1e-8
		Y_perturb = np.copy(Y_guess)
		F_perturb = np.zeros((3,1))
		Jac = np.zeros((3,3))
		k=0
		for i in range(n):

			Y_perturb[i] = Y_guess[i] + eps
			F_perturb[0] = f1(Y_old[0], Y_perturb[0], Y_perturb[1], Y_perturb[2], dt)
			F_perturb[1] = f2(Y_old[1], Y_perturb[0], Y_perturb[1], Y_perturb[2], dt)
			F_perturb[2] = f3(Y_old[2], Y_perturb[0], Y_perturb[1], Y_perturb[2], dt)
			Jac[0,i] = (F_perturb[0] - F[0])/eps
			Jac[1,i] = (F_perturb[1] - F[1])/eps
			Jac[2,i] = (F_perturb[2] - F[2])/eps
			#print(X_perturb)
			Y_perturb[i] = Y_guess [i]
		#print(Jac)

		#Computing the new values
		Y_new = Y_guess - alpha*np.matmul(inv(Jac),F)

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

		#Updating old values 
		Y_guess = Y_new

		#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)
		iter = iter + 1

	#Final results
	y1 = Y_new[0]
	y2 = Y_new[1]
	y3 = Y_new[2]

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

	time_iter = time_iter + 1 

#Make sure the size of the matrix and time vector corresponds
size = np.shape(Y)
print(size)
print(total_iter)

plt.plot(t,Y[0,:],color='black')
plt.plot(t,Y[1,:],color='orange')
plt.plot(t,Y[2,:],color='blue')
plt.legend(['Y1','Y2','Y3'])
plt.show()

Projects by Jorge Martinez

The purpose of this project is to show the effect of initial temperature and pressure in the ignition delay in a combustion process. For this purpose, a Constant Volume Reactor and the corresponding network objects are created in Cantera, and the combustion process is s Read more

The purpose of this project, is to setup and analyze a transient shock tube simulation using Converge Studio. The project will analyze how the interaction of the original shock waves and the reflected waves affect the fluid in the tube, and how the computational mesh ha Read more

The purpose of this project is to simulate the conjugate heat transfer in a pipe using Converge Studio. For this, an air flow is simulated through an aluminum pipe, which is receiveing a constant heat flux of 10,000 W/m2 at the outside wall. The objective is to understa Read more

This project objective is to simulate the Prandtl Meyer shock wave phenomena using Converge Studio. First, a quick literature review about shock waves and their boundary conditions is provided. Then , the problem is set up and solved, focusing in how the different param Read more

The objective of this part of the project, is to simulate the transient flow through an elbow pipe with a throttle valve in the middle. This valve will be rotating from zero to 25 degrees, and then stay steady.  Case Setup. In the steady state part of the project Read more

This project purpose is to simulate the flow over an elbow body that contains a throttle using CONVERGE Studio. In this first part of the project, an steady state analysis will be set up. This means that the throttle will not move. This case helps set up the real analys Read more

In this challenge, the efficiency of using a heat recuperator to pre-heat the air is analyzed. For that, we simulate a furnace burning methane and air mixture, and we analyze first how the pre-heating affects at the AFL, and then the total heat/work that could be extrac Read more

This project main objective is to calculate the Adiabatic Flame Temperature using Python and Cantera. In the first part of the project, the effect of the equivalence ratio for the methane combustion in a constant volume chamber is analyzed, and the python and cantera re Read more

This projects aims to simulate a backward facing step flow using CONVERGE Studio. The flow is creating when a not reacting specie, air (23% O2 and 77% N2 in mass fraction), suffers a pressure gradient between the ends of the channel. The velocity, pressure and recirulat Read more


Loading...

The End