Week 6 - Multivariate Newton Rhapson Solver

Objective: To solve a given set of Ordinary Differential equations using the Multi-Variate Newton Raphson Method.

Given: The set of ODE's are given below:

`dy_1/dt= - 0.04*y_1+10^4*y2*y3`

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

The jacobian should be estimated numerically and not analytically. The initial conditions are 1 for y1 and 0 for y2 and y3 respectively.
Take alpha as 1. Run the simulation for 10 minutes.

Theory:

                                                                                                           `"Multi-Variate Newton Raphson Method"`

 As we know, the Multivariate Newton-Raphson method is a direct extension of the single variable Newton-Raphson method. Where the single variable Newton-Raphson method solved
`f(x)=0`, the multivariate version will solve a system of n equations of the form:

`f_1(x_1,x_2,x_3.....,x_(n-1),x_n) = 0`

`f_2(x_1,x_2,x_3.....,x_(n-1),x_n) = 0`

`f_3(x_1,x_2,x_3.....,x_(n-1),x_n) = 0`
.
.
.
.

`f_(n-1)(x_1,x_2,x_3.....,x_(n-1),x_n) = 0`

`f_n(x_1,x_2,x_3.....,x_(n-1),x_n) = 0`

We will adopt the shorthand notation for the equation.

`f(x)=0`

According to Multi-variate Newton Raphson Method, to get the solution out of coupled non-linear ODE below equation we use:

`y_("new")=y_("old")-alpha*J^(-1)*F`

where, alpha - Relaxation factor

Vector F will be written as,

`f=[[f_1],[f_2],[f_3],[.],[.],[f_n]]`

Old value or Guess value will be written as,

`y_("old")=[[y_("1 old")],[y_("2 old")],[y_("3 old")],[.],[.],[y_("n old")]]`

Jacobian Matrix,

`J=[[(df_1)/dy_1(df_1)/dy_2(df_1)/dy_3],[(df_2)/dy_1 (df_2)/dy_2(df_2)/dy_3],[(df_3)/dy_1(df_3)/dy_2(df_2)/dy_3],[.],[.]]`

where, nos. of column= Nos. of independent variables
           nos. of Rows= Nos. of Functions

Solution:

Python Code:

"""
Program to solve non-linear coupled ODE with multi-variate Newton Raphson Solver

System of coupled non-linear ODE is:

  -0.04 * Y1 + 10^4 * Y2*Y3 = 0
  0.04 * Y1 - 10^4 *Y2*Y3 - 3 * 10^7* Y2^2 = 0
  3 * 10^7* Y2^2 = 0

By- Gaurav V Kharwade || Skill-Lync:2020
"""

import numpy.matlib 
import numpy as np
import numpy.linalg
import matplotlib.pyplot as plt

def f1(y1,y2,y3,y1_old,dt): # First non-linear equation
	return (y1_old + dt*(-0.04*y1 + 1e4*y2*y3) - y1)

def f2(y1,y2,y3,y2_old,dt): # Second non-linear equation
	return (y2_old + dt*(0.04*y1 - (1e4*y2*y3) - (3*1e7*y2*y2)) - y2)

def f3(y1,y2,y3,y3_old,dt): # Third non-linear equation
	return (y3_old + dt*(3e7*y2*y2) - y3)

def jacobian(y1,y2,y3,y1_old,y2_old,y3_old,dt): # To define Jacobian matrix

	#To define 3 * 3 matrix creation
	j= np.ones((3,3))
	
	# Step-size
	h = 1e-6
	
	"""
	For estimation of jacobian matrix we have used numerical method(BDS) instead analytical method
	"""
	# Row 1
	j[0,0]= (f1(y1+h,y2,y3,y1_old,dt)-f1(y1,y2,y3,y1_old,dt))/h
	j[0,1]= (f1(y1,y2+h,y3,y1_old,dt)-f1(y1,y2,y3,y1_old,dt))/h
	j[0,2]= (f1(y1,y2,y3+h,y1_old,dt)-f1(y1,y2,y3,y1_old,dt))/h

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

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

	return j

# Initialisation
y1_old= 1
y2_old= 0
y3_old= 0

# Old value assignment
y_old= np.ones((3,1))
y_old[0]= y1_old
y_old[1]= y2_old
y_old[2]= y3_old

# Column vector of function F
F= np.copy(y_old)

# Relaxation factor
alpha= 1
# Total time, T=10 mins, t= 10 * 60 
t= 600
# Time step
dt= 0.1
n= int(t/dt)

# Array of updated values
y_new1= []
y_new2= []
y_new3= []
it= []
error_updated= []

# Initialising iteration
iter= 1

# For loop, to time stepping
for i in range(0,n):
	
	# Newton Raphson criteria
	tol= 9e-6
	error= 9e9

	# Newton Raphson Solver
	while error > tol:

		# Jacobi Matrix
		j= jacobian(y_old[0],y_old[1],y_old[2],y1_old,y2_old,y3_old,dt)

		# F column vector
		F[0]= f1(y_old[0],y_old[1],y_old[2],y1_old,dt);
		F[1]= f2(y_old[0],y_old[1],y_old[2],y2_old,dt);
		F[2]= f3(y_old[0],y_old[1],y_old[2],y3_old,dt);

		# Equation for solver
		y_new = y_old - alpha * np.matmul(numpy.linalg.inv(j), F);

		# Error finding
		error = np.max(np.abs(y_new - y_old));
		print(error)

		# Old value updation
		y_old = y_new
			
	# To print on console
	log= 'iter= {0} error= {1} y1= {2} y2= {3} y3= {4} '.format(iter,error, y_new[0], y_new[1], y_new[2])
	print(log)

	y_old = y_new;
	y_new1.append(y_new[0])
	y_new2.append(y_new[1])
	y_new3.append(y_new[2])
	error_updated.append(error)

	
	iter= iter + 1
	it.append(iter)	
	
	# Old value updation
	y1_old= y_new[0]
	y2_old= y_new[1]
	y3_old= y_new[2]

# To create a plot of error dropping	
plt.plot(it,error_updated,'--',color='red')
plt.yscale("log")
plt.title('Error dropping characteristics')
plt.xlabel('Iteration nos.')
plt.ylabel('Error value')
plt.grid('on')

# To create a figure with 3 rowa and 1 column of subplot
fig, ax= plt.subplots(3)
# For convergence of y1
ax[0].plot(it,y_new1)
ax[0].set_title('Convergence of y1')
ax[0].grid('on')
# For convergence of y1
ax[1].semilogy(it,y_new2)
ax[1].set_title('Convergence of y2')
ax[1].grid('on')
# For convergence of y3
ax[2].plot(it,y_new3)
ax[2].set_title('Convergence of y3')
ax[2].grid('on')

fig.tight_layout()
plt.show()

Results:

After running the program the results we get is as below:

We can see above the error values that we have targeted have achieved and the solution for this ODE has shown above.
y1= 0.39997579, y2= 2.63177e-6 & y3= 0.600021

We can see below, at each and every iterations we are getting closer to an accurate value.

In this plot, we can see how the errors are getting dropped at every iteration in order to achieve a solution with high accuracy.

I ran the program by using a time step much lower than 0.1, but it has been observed that there was no major variation in the solution obtained. Hence, choosing 0.1 as time step will be a better choice as it helps us to get solutions in less time and with high accuracy also.