Week 6 - Multivariate Newton Rhapson Solver

In this challenge, you will solve the set of equations given below,

`dy_1/dt = -0.04y_1 + 10^4y_2y_3` .......(1)

`dy_2/dt = 0.04y_1-10^4y_2y_3-3.10^7y_2^2` ........(2)

`dy_3/dt = 3 . 10^7y_2^2` .......(3)

Solve the problem using Implicit Euler Method/Backward Differencing, assume the initial conditions as 1 for y1 and 0 for y2 and y3 respectively.

Take alpha as 1. Run the simulation for 10 minutes.

solution:

Theory :

After applying the Backward Differencing Method to the system of ODE, the equation becomes a coupled non-linear differential equation and as there is more equation, Multivariate Newton Raphson method is used to solve the above
equation

Newton-Raphson Method (Multivariate) (*http://fourier.eng.hmc.edu/e176/lectures/ch2/node6.html)

The Newton-Raphson method discussed above for solving a single-variable equation can be generalized for solving multivariate equation systems containing equations of variables in :

we first consider the Taylor series expansion of each of the functions of variables in around a point :

where represents the second and higher order terms in the series beyond the linear term, which can be neglected if is small. These equations can be expressed in matrix form

`f(x_0) + J(x_0) (x-x_0) +r approx f(x_0) + J(x_0) (x-x_0)`

where is an Jacobian matrix of evaluated at .

for given equation system f(x) = 0 , we have M = N : The number of equations is the same as the number of unknowns, the Jacobian J (x) is a square matrix and in general invertible with an inverse `J^-1` exists. In the special case where all functions in f(x) are linear, the Taylor expansion above contains only the first two terms (constant and linear) while all higher order terms are zero. The approximation becomes exact. To find the root at which , we set the above function to zero f(x) = 0, and solve the resulting equation for x to get.

`x^* = x = x_0 - J(x_0)^-1 f(x_0)`

which can be improved iteratiively to move from an initial point towards the root :

Same as in the univariate case, a scaling factor `alpha` can be used to control the step size of the iteration

`x_(n+1) = x_n + alpha Deltax_n = x_n - alpha J_n^-1 f_n`

When ,`alpha < 1` the step size becomes smaller and the convergence of the iteration is slower, however, we will have a better chance not to skip a solution, which may happen if f (x) is not smooth and the step size is too big.

The algorithm is listed below:

• Select `X_0`
• Obtain`f_0 = f(x_0) and J_0`
• Obtain `J_0^-1`
• n = 0
• While do
  • `x_(n+1) = x_n - J_n^-1 f_n`
  • find `f_(n+1) and J(n+1)`
  • n = n+1

Backward Differencing Method:

`dy/dt = (y^n -y^(n-1))/(Deltat)`

so above set of equations we can write as given below,

`(y_1^n -y_1^(n-1))/(Deltat) = -0.04y_1 + 10^4y_2y_3``

`(y_2^n -y_2^(n-1))/(Deltat)= 0.04y_1-10^4y_2y_3-3.10^7y_2^2`

`(y_3^n -y_3^(n-1))/(Deltat)= 3 . 10^7y_2^2`

`y_1^(n-1) + Deltat (-0.04y_1^n + 10^4y_2^ny_3^n) - y_1^n = 0`

`y_2^(n-1) + Deltat ( 0.04y_1-10^4y_2y_3-3.10^7y_2^2) - y_2^n = 0`

`y_3^(n-1) + Deltat (3 . 10^7y_2^2) - y_3^n = 0`

in above equations `y^n is y^(New)` & `y^(n-1) is y^(old)`

Python code :

"""
Program to solve the system of coupled non-linear equations
dy1/dt = -0.04*y1 + 10^4*y2*y3
dy2/dt = 0.04*y1 - 10^4*y2*y3 - 3*10^7*y2*y2
dy3/dt = 3*10^7*y2*y2
Using Newton Raphson Multivariate Method
"""
import numpy as np
from numpy.linalg import inv
import matplotlib.pyplot as plt

def f1(y1_old, y1, y2, y3, dt):			#1st non-linear equation
	
	return (y1_old + dt*(-0.04*y1 + y2*y3*pow(10,4)) - y1)

def f2(y2_old, y1, y2, y3, dt):			#2nd non-linear equation
	
	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):			#3rd non-linear equation

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

def jacobian(y1_old, y2_old, y3_old, y1, y2, y3, dt):		#Defining the Jacobian Matrix

	J = np.ones((3,3))		# 3 x 3 unit matrix

	h = 1e-8		#Step Size
	#Row 1
	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
	#Row 2
	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
	#Row 3
	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

# initial conditions values

y1_old = 1
y2_old = 0
y3_old = 0

# Initial error

error = 9e9

#numpy column vector

Y_old = np.ones((3,1))

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

#numpy column vector for F
F = np.copy(Y_old)

#parameters for newton raphson calculations

alpha = 1
tol = 1e-9
t = 10*60
dt = 0.5
n = int(t/dt)

y_1 = []
y_2 = []
y_3 = []

iteration = []
error_count = []
sub_iteration = []
#starting newton raphson loop
sub_iter = 1
iter = 1

for i in range (0,n):
	while error > tol:
		# Computing the jacobian old values are used
		J = jacobian(y1_old,y2_old,y3_old, Y_old[0], Y_old[1], Y_old[2], dt)

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

		#Computing new values
		Y_new = Y_old - alpha*np.matmul(inv(J), F)

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

		#Updating old values
		Y_old = Y_new
		sub_iter = sub_iter + 1



	#log message
	log_message = 'iteration number = {0} sub_iter =  {1} y1 = {2} y2 = {3} y3 =  {4} error = {5} dt {6}'.format(iter, sub_iter,Y_new[0], Y_new[1], Y_new[2] , error ,dt)
	print(log_message)
	iter = iter+1

	#Storing values for plotting
	y_1.append(Y_new[0])
	y_2.append(Y_new[1])
	y_3.append(Y_new[2])
	iteration.append(iter)
	sub_iteration.append(sub_iter)
	error_count.append(error)
	Y_old = Y_new
	y1_old = Y_new[0]
	y2_old = Y_new[1]
	y3_old = Y_new[2]
	error = 9e9

#Final Results

plt.plot(iteration,y_1)
plt.plot(iteration,y_2)
plt.plot(iteration,y_3)
plt.plot(iteration,error_count)
#plt.plot(iteration,sub_iteration)
plt.xlabel('Number of steps')
plt.ylabel('Variables(Species)')
plt.grid ('on')
plt.legend(['y_1', 'y_2', 'y_3','error_count'])
#plt.legend(['error_count'])
#plt.legend(['sub_iteration'])
plt.title('Multivariate Newton Raphson Method for dt = 0.5')
plt.show()

Choosing a suitable time step size:-

The Python code has been run for different time step that 1, 0.1,0.5 & 0.01 receptively. But there is no significant change in the values of y1,y2,y3. As we have implemented implicit method , no significant effect of time-step on results. Total time for calculations are 4.4s ,10.1s,7.4s , 55.09s. So we may have take time-step as 0.5 or 0.1 as they have good range of error of order e-10 and e-12 or we can have variable time-step if we have complex equation.