Multivariate Newton Rahpson method   (Cantera, Skill Lync, Medium)

Aim:

To understand and create a program for Multivariate Newton Raphson method.

Introduction:

In many areas of engineering, especially in combustion simulations, the equations which are to be solved are differential equations. Even when the problem under consideration will be small, the equations will be interdependent or coupled. Solving these equations analytically is not possible, so they are numerically integrated using Newton Raphson method. It is a multivariate method, the univariate method is its derived form for single variable. In this report, we solve a set of coupled equations using multivariate method and implicit formulation.

The differential equations which are solved are:

The initial conditions are: y1 = 1, y2 = 0, y3 = 0

The time-differential terms are discretized using the euler's backward difference scheme or implicit method.

The superscript denotes the time-step at which the values are taken.

The jacobian matrix is defined as:

In this report, the differential terms of jacobian matrix are calculated using numerical approximation:

The F vector is the matrix of function (equations), which are to solved for.

The values for next timestep are found using this formula, which is the original formula of multivariate newton raphson method.

How does the code work?

1) The values of time for which solution is to be calculated is are stored in an array.

2) The variables for initial solution, old values, new values and guess values are initialized.

3) The initial error and tolerance values are stated

4) Outer loop begins, which calculates the solution for each time.

5) The inner loop begins after few updates in variable data. This inner loop uses the multivariate newton raphson method. It continues till error is larger than tolerance

6) The jacobian matrix is calculated using function 'jacob'

7) The values of F vector are calculated using f1, f2 and f3 functions.

8) The newton raphson formula is used and required updates in variable values are done.

9) After the loops end, the solution is plotted.

"""
Program to solve a system of ODEs using Multivariate Newton-Raphson method

current example ODEs:
	d(y1)/dt = -0.04*y1 + (10^4)*y2*y3
	d(y2)/dt = 0.04*y1 - (10^4)*y2*y3 - 3*(10^7)*(y2)^2
	d(y3)/dt = 3*(10^7)*(y2)^2

"""

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

#Functions that would be used to find roots
def f1(y1_o,y1_n,y2_n,y3_n,dt):
	return y1_o - 0.04*dt*y1_n + 1e4*dt*y2_n*y3_n - y1_n

def f2(y2_o,y1_n,y2_n,y3_n,dt):
	return y2_o + 0.04*dt*y1_n - 1e4*dt*y2_n*y3_n - 3e7*dt*pow(y2_n,2) - y2_n

def f3(y3_o,y1_n,y2_n,y3_n,dt):
	return y3_o + 3e7*dt*pow(y2_n,2) - y3_n

#Jacobian matrix
def jacob(nY,oY,dt):
	dy = 1e-5
	J = np.zeros((3,3))

	J[0,0] = (f1(oY[0],nY[0]+dy,nY[1],nY[2],dt) - f1(oY[0],nY[0],nY[1],nY[2],dt))/dy
	J[0,1] = (f1(oY[0],nY[0],nY[1]+dy,nY[2],dt) - f1(oY[0],nY[0],nY[1],nY[2],dt))/dy
	J[0,2] = (f1(oY[0],nY[0],nY[1],nY[2]+dy,dt) - f1(oY[0],nY[0],nY[1],nY[2],dt))/dy

	J[1,0] = (f2(oY[1],nY[0]+dy,nY[1],nY[2],dt) - f2(oY[1],nY[0],nY[1],nY[2],dt))/dy
	J[1,1] = (f2(oY[1],nY[0],nY[1]+dy,nY[2],dt) - f2(oY[1],nY[0],nY[1],nY[2],dt))/dy
	J[1,2] = (f2(oY[1],nY[0],nY[1],nY[2]+dy,dt) - f2(oY[1],nY[0],nY[1],nY[2],dt))/dy

	J[2,0] = (f3(oY[2],nY[0]+dy,nY[1],nY[2],dt) - f3(oY[2],nY[0],nY[1],nY[2],dt))/dy
	J[2,1] = (f3(oY[2],nY[0],nY[1]+dy,nY[2],dt) - f3(oY[2],nY[0],nY[1],nY[2],dt))/dy
	J[2,2] = (f3(oY[2],nY[0],nY[1],nY[2]+dy,dt) - f3(oY[2],nY[0],nY[1],nY[2],dt))/dy

	return J

#Constructing t array
t_start = 0
t_stop = 600
dt = 1
n = (t_stop-t_start)/dt +1
t = np.linspace(t_start, t_stop,int(n))

oY = np.zeros((3,1))
oY[0] = 1
F = np.copy(oY) #Using Deep copy method instead of passing address.
nY = np.zeros((3,1))

err = 9e9
alpha = 1
tol = 1e-7
count = 0
y1 = [1]*len(t)
y2 = [0]*len(t)
y3 = [0]*len(t)

for i in range(1,len(t)):
	y1_o = oY[0]
	y2_o = oY[1]
	y3_o = oY[2]
	gY = np.copy(oY)

	while err>=tol:

		jac = jacob(nY,oY,dt)
		F[0] = f1(y1_o,gY[0],gY[1],gY[2],dt)
		F[1] = f2(y2_o,gY[0],gY[1],gY[2],dt)
		F[2] = f3(y3_o,gY[0],gY[1],gY[2],dt)

		nY = gY - alpha*np.matmul(inv(jac),F)
		err = max(abs(nY - gY))

		gY = nY
		count += 1

	oY = nY
	
	y1[i] = nY[0]
	y2[i] = nY[1]
	y3[i] = nY[2]


	log_message = 'Time = {0} sec, iterations = {1}, y1 = {2}, y2 = {3}, y3 = {4}'.format(round(t[i],5),count,nY[0],nY[1],nY[2])
	print(log_message)
	count = 0
	err = 9e9


plt.plot(t,y1,label = '$y1$')
plt.plot(t,y2,label = '$y2$')
plt.plot(t,y3,label = '$y3$')
plt.xlabel('Time (sec)')
plt.ylabel('$y_i$')
plt.title('Multivariate Newton-Raphson Solution for dt = '+ str(dt)+ ' sec')
plt.legend()
plt.grid(True)
plt.show()

Results:

As it can be seen, all the timesteps approximately give the same values. This is due to use of implicit formulation. It allows us to choose higher time step sizes. The ideal time step size depends on how much accurate results we want. For this we can solve the system for some range of timesteps and find the difference in their values. The timestep which would give the required accuracy will be the ideal time step.

Conclusion:

The method of Multivariate Newton Raphson and use of implicit method is understood.