Week 6 - Multivariate Newton Rhapson Solver

Objectve : To solve the given coupled non linear system of equations using mulltivariate newton raphson method.

Theory

The given set of differential equations are as follows:

The above differential equations using backward differencing scheme can be written as:

Writing the above equations in the form of root finding problem so that multivariate Newton Raphson method can be used.

The initial guess values for y1, y2 and y3 are taken as 1, 0 and 0 respectively.

Jacobian matrix

Jacobian matrix is the collection of all the possible derivatives for a given system of equations

Then mulltivariate newton raphson method is used to compute new values till error value is greater than the given tolerance.

Y^new = Y^old – α.J^-1.F

Where,

Y^new is the new iteration value calculated from initial guess value

Y^old is the value obtained from the previous iteration.

α is the relaxation factor

F is the matrix containing the given system of equations.

J^-1 is the inverse of the Jacobian matrix

Python program

"""
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

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

def f1(y1_old, y1, y2, y3, dt):

	return y1_old + dt*(-0.04*y1 + (pow(10,4)*y2*y3)) - y1

def f2(y2_old, y1, y2, y3, dt):

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

def f3(y3_old, y1, y2, y3, dt):

	return y3_old + dt*(3*pow(10,7)*pow(y2,2)) - y3

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

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

	 h = 1e-5

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

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

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

	 return J

# Defining the initial guess values

y1_old = 1
y2_old = 0
y3_old = 0

error = 9e9

Y_old = np.ones((3,1))
Y_old[0] = y1_old
Y_old[1] = y2_old
Y_old[2] = y3_old

F = np.copy(Y_old)

alpha = 1
tol = 1e-6
t = 600
dt = 0.1
n = int(t/dt)
y_1 = []
y_2 = []
y_3 = []
iteration = []

iter = 1


for i in range (0, n):

	while error > tol:

		J = jacobian(y1_old, y2_old, y3_old, Y_old[0], Y_old[1], Y_old[2], dt)

		F[0] = f1(y1_old, Y_old[0], Y_old[1], Y_old[2], dt)
		F[1] = f2(y2_old, Y_old[0], Y_old[1], Y_old[2], dt)
		F[2] = f3(y3_old, 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)

		Y_old = Y_new

	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

	#plotting

	y_1.append(Y_new[0])
	y_2.append(Y_new[1])
	y_3.append(Y_new[2])

	iteration.append(iter)

	Y_old = Y_new
	y1_old = Y_new[0]
	y2_old = Y_new[1]
	y3_old = Y_new[2]

	error = 9e9

plt.plot(iteration, y_1, color = 'r')
plt.plot(iteration, y_2, color = 'g')
plt.plot(iteration, y_3, color = 'b')
plt.xlabel ('Time (s)')
plt.ylabel ('Species (Y)')
plt.grid('on')
plt.legend(['y1', 'y2', 'y3'])
plt.title ('Multivariate Newton Raphson Solver')
plt.show()

Plots for different cases

Case 1 : dt = 0.01 (n = 60000)

Case 2 : dt = 0.1 (n = 6000)

Case 3 : dt = 1 (n = 600)

Conclusion

For the time step dt, we must consider a value which should be able to give accurate results but at the same time it should not take too much computation time so that higher computation power is not required. When we take dt as 1, there are 600 steps which is large for this case and hence the results are not accurate. When take we take dt as 0.01, there are 60000 time steps. In this case the results are accurate but the time taken for computation is long. Hence for this case, the preferred value of dt is 0.1 which has 6000 time steps and gives results which are accurate with less computation time.