Week 6 - Multivariate Newton Rhapson Solver

Aim :To solve the set of equations with multiple variables using the Newton Raphson Solver.

Apparatus or given data:

1)`(d underset(1)(y))/dt=-0.04underset(1)(y)+10^4underset(2)(y)underset(3)(y)`

2)`(d underset(2)(y))/dt=0.04underset(1)(y)-10^4underset(2)(y)underset(3)(y)-3.10^7underset(2)(y^2)`

3)`(d underset(3)(y))/dt=3.10^7underset(2)(y^2)`

Initial Conditions:

`underset(1)(y)=1; underset(2)(y)=0; underset(3)(y)=0;`

Methods for Numerical analysis:

1)Implicit Euler method

2)Jacobian method

Solver :

Newton Raphson solver

Procedure:

The set of equations are

`(d underset(1)(y))/dt=-0.04underset(1)(y)+10^4underset(2)(y)underset(3)(y)`

`(d underset(2)(y))/dt=0.04underset(1)(y)-10^4underset(2)(y)underset(3)(y)-3.10^7underset(2)(y^2)`

`(d underset(3)(y))/dt=3.10^7underset(2)(y^2)`

The change here is that we are incorporating a new step process to solve the multivariate system of equations as it is not possible to solve these equations using the normal newton raphson technique.

The Newton-Raphson method (also known as Newton's method) is a way to quickly find a good approximation for the root of a real-valued function f(x) = 0. It uses the idea that a continuous and differentiable function can be approximated by a straight line tangent to it.

Suppose you need to find the root of a continuous, differentiable function f(x), and you know the root you are looking for is near the point . Then Newton's method tells us that a better approximation for the root is

`underset(1)(x)=underset(0)(x)-[(f(underset(0)(x)))/(f'(underset(0)(x)))].a`

where, a = relaxation factor

This process may be repeated as many times as necessary to get the desired accuracy. In general, for any x-value, the next value is given by,

`underset(n+1)(x)=underset(n)(x)-[(f(underset(n)(x)))/(f'(underset(n)(x)))].a`

For this, a guess value, x0 is first taken into consideration using a ball park number which is approximately closer to the roots for the given function. Then by selecting an optimum relaxation factor, we can determine the accurate roots of the function without too many iterations than necessary. This same method used for the univariate method is developed for the multivariate solver as,

`underset(New)(x)=underset(guess)(x)-alpha.j^-1.F`

where,

J = Jacobian matrix of the system of equations.

The Jacobian matrix is basically a matrix consisting the first order partial derivatives of the multivariate system of equations.

J=`[[((dunderset(1)(F))/(dunderset(1)(y))) ((d underset(1)(F))/(dunderset(2)(y))) ((dunderset(1)(F))/(dunderset(3)(y)))],[((dunderset(2)(F))/(dunderset(1)(y))) ((dunderset(2)(F))/(dunderset(2)(y))) ((dunderset(2)(F))/(dunderset(3)(Y)))],[((dunderset(3)(F))/(dunderset(1)(y))) ((dunderset(3)(F))/(dunderset(2)(y))) ((dunderset(3)(F))/(dunderset(3)(y)))]]`

where,

F is the function matrix, F = `[[f1],[f2],[f3]]`

DERIVATION

We know that the system of equations are,

`(d underset(1)(y))/dt=-0.04underset(1)(y)+10^4underset(2)(y)underset(3)(y)` ;

`(d underset(2)(y))/dt=0.04underset(1)(y)-10^4underset(2)(y)underset(3)(y)-3.10^7underset(2)(y^2`

`(d underset(3)(y))/dt=3.10^7underset(2)(y^2)`

The initial conditions given are,

y1 = 1

y2 = 0

y3 = 0

Now, all the equations are discretized using the backward differencing as,

`(underset(1)(y^n)-underset(1)(y^(n-1)))/(gradt) = -0.04underset(1)(y^n)+10^4underset(2)(y^n)underset(3)(y^n)`

`(underset(2)(y^n)-underset(2)(y^(n-1)))/(gradt) = 0.04underset(1)(y^n)-10^4underset(2)(y^n)underset(3)(y^n)-3.10^7(underset(2)(y^n))^2`

`(underset(3)(y^n)-underset(3)(y^(n-1)))/(gradt) = 3.10^7(underset(2)(y^n))^2`

Simplified further, we get a set of equations,

`underset(1)(y^n)=underset(1)(y^(n-1))+h.(-0.04underset(1)(y^n)+10^4underset(2)(y^n)underset(3)(y^n))`

`underset(2)(y^n)=underset(2)(y^(n-1))+h.(0.04underset(1)(y^n)-10^4underset(2)(y^n)underset(3)(y^n)-3.10^7(underset(2)(y^n))^2)`

`underset(3)(y^n)=underset(3)(y^(n-1))+h.(3.10^7(underset(2)(y^n))^2)`

where, h = dt

After this, we write the equation in root finding function form, i.e, we can equate the terms to zero, and hence get,

`underset(1)(f)=underset(1)(y^(n-1))+h.(-0.04underset(1)(y^n)-10^4underset(2)(y^n)underset(3)(y^n)-underset(1)(y^n))=0`

`underset(2)(f)=underset(2)(y^(n-1))+h.(0.04underset(1)(y^n)-10^4underset(2)(y^n)underset(3)(y^n)-3.10^7(underset(2)(y^n))^2)-underset(2)(y^n)=0`

`underset(3)(f)=underset(3)(y^(n-1))+h.(3.10^7(underset(2)(y^n))^2)-underset(2)(y^n)=0`

Hence, using this function matrix, we will calculate the Jacobian matrix and then solve using the Multivariate Newton raphson solver

`underset(New)(Y)=underset(guess)(Y)-alpha.(j^(-1).F)^n`

Coding:

For solving the actual problem using Python, we proceed as usual and solve the quadratic equations first by defining functions for each equation. The Jacobian matrix is then differentiated numerically, by discretizing using the forward differencing scheme. An arbitraty dy value is assumed for the accurate discretization of the PDE. In this case, dy = 1e-10.

Dummy values are assumed in an array which will be used in the time loop required for the problem. In our case, we are asked to simulate the system for 10 mins, i.e, for 600 seconds. So the solution arrays defined are given lengths according to the total time and the time step required to solve the system.

Python Code:

'''
Program to solve a multivariate system using Newton Raphson solver
'''
import numpy as np 
import matplotlib.pyplot as plt
from numpy.linalg import inv

# Functions to solve each multivariable equation
def f1(y1,y2,y3,y_old1,h): 
	return y_old1 + h*(-0.04*y1 + (pow(10,4)*y2*y3)) - y1

def f2(y1,y2,y3,y_old2,h): 
	return y_old2 + h*(0.04*y1 - (pow(10,4)*y2*y3) - (3*pow(10,7)*pow(y2,2))) - y2

def f3(y1,y2,y3,y_old3,h): 
	return y_old3 + h*(3*pow(10,7)*pow(y2,2)) - y3

# function to formulate the jacobian matrix for the system of equations
def jacobian(y_new,y_old,h):

	# input an identity 3x3 matrix
	J = np.ones((3,3))

	# Using dummy values that will be updated for each time step when solving
	y1 = y_new[0] 
	y2 = y_new[1] 
	y3 = y_new[2] 

	# Old values are the input guess values we give
	y_old1 = y_old[0]
	y_old2 = y_old[1]
	y_old3 = y_old[2]

	# To solve the system numerically
	dy = 1e-10

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

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

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

	return J

# Input the guess values as old, due to use of backward differencing
y_old = np.zeros((3,1))
y_old[0] = 1

y_new = np.zeros((3,1))

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

# inputs
a = 1 # alpha for Newton raphson solver
tol = 1e-10
error = 1e9

dt = 10 # time step
h = dt

# time to run simulation = 10 mins = 600 secs
t = np.arange(0,600,h)

# Empty arrays initialized for storing final values
e = np.zeros((int(600/h),1))
Y1 = np.zeros((int(600/h),1))
Y2 = np.zeros((int(600/h),1))
Y3 = np.zeros((int(600/h),1))

# time loop started
for i in range(0,len(t)):
	y_old1 = y_old[0]
	y_old2 = y_old[1]
	y_old3 = y_old[2]
	y_guess = np.copy(y_old)

	while error > tol:
		# computing jacobian - old values used
		J = jacobian(y_new,y_old,h)

		y1g = y_guess[0]
		y2g = y_guess[1]
		y3g = y_guess[2]

		# Computing F vector - old values used
		F[0] = f1(y1g,y2g,y3g,y_old1,h)
		F[1] = f2(y1g,y2g,y3g,y_old2,h)
		F[2] = f3(y1g,y2g,y3g,y_old3,h)

		# computing new values
		y_new = y_guess - a*np.matmul(inv(J),F)
		error = max(abs(y_new-y_guess))

		y_guess = y_new
	
	# Values after each time step added into array
	e[i] = error
	y_old = y_new
	Y1[i] = y_new[0]
	Y2[i] = y_new[1]
	Y3[i] = y_new[2]
	error = 1e9

# Plots
plt.plot(t,Y1, label = 'y1')
plt.plot(t,Y2, label = 'y2')
plt.plot(t,Y3, label = 'y3')
plt.title('Multi-variate Newton Raphson Solver for dt = ' +str(h))
plt.xlabel('Time (secs)')
plt.ylabel('$y_i$')
plt.legend()
plt.grid()
plt.show()

Results:

About the code :

The error and initial values are defined before the time loop starts and are accordingly updated at the end of the loop for the next iteration to take place. The 'matmul' function from Python is used to facilitate the matrix multiplication needed for the solver.

Result Explanation:

For the time step dt, we must assume the best value that gives a fairly accurate result for the problem. If we take dt as 1, there are 600 steps in the system and hence results are sharp and not that accurate. But if we give a smaller dt of say, 1e-2, we have too many iterations to complete in an already complex multivariate solver, eventhough the solution result comes out smooth. Hence, an average dt value of 0.1 is assumed in this case to give an accurate representation of data and also sufficiently time efficient.

The results show the variations in values and their domains getting smaller due to increase in time step size. Hence the ideal time step for this case is assumed to be 0.1.