Week 6 - Multivariate Newton Rhapson Solver

Aim

To solve a set of coupled ordinary differential equations using multivariate Newton Rhapson solver

1. Introduction

For an ordinary differential equation with independent variable t and dependent variable y, using Newton Rhapson method, the new value can be calculated as

y_new = y_guess – `f(y_(guess),t)/f^'(y_(guess),t)`

For a coupled ordinary differential equations with multiple variables, we need to modify this solver by incorporating column matrices for the guess values and function values. Furthermore Jacobian matrix is included to arrive at the equation

y_new = y_guess – `alpha*J^(-1).F`

Where,

`alpha` = relaxation factor

J= Jacobian matrix

F= Function value at new step size

For a case of ODEs with three dependentvariables y1, y2 and y3 defined with independent variable t, we can write,

Ynew = `[(y_1^(n+1)),(y_2^(n+1)),(y_3^(n+1))]`

Yold= `[(y_1^(n)),(y_2^(n)),(y_3^(n))]`

F=`[(f_1^(n+1)),(f_2^(n+1)),(f_3^(n+1))]`

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_3)/dy_3]]`

2. Objective

To solve the below set of equations using multivariate NR method and incorporate implicit Euler/backward differencing method to arrive at a stable solution (and address the stiffness of solution)

`dy_1/dt = -0.04 y_1 +10^4.y_2.y_3`

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

`dy_3/dt = 3.10^7.y_2^2`

3. Equations

Using implicit Euler method, we can rewrite the ordinary differential equations as

`y_1^n = y_1^(n-1) + dt(-0.04 y_1 +10^4.y_2.y_3)^n`

`y_2^n = y2^(n-1) + dt(0.04y_1 - 10^4.y_2.y_3 - 3.10^7y_2^2)^n`

`y_3^n =y_3^(n-1) + dt(3.10^7.y_2^2)`

By subtracting RHS-LHS, we get

`y_1^(n-1) + dt(-0.04 y_1 +10^4.y_2.y_3)^n -y_1^n = 0 = f_1`

`y2^(n-1) + dt(0.04y_1 - 10^4.y_2.y_3 - 3.10^7y_2^2)^n-y_2^n=0=f_2`

`3.10^7.y_2^2-y_3^n = 0= f3`

4. Program algorithm

The following are the steps followed in the program:

1. To calculate the matrix operations, we import numpy library and from numpy linear algerbra, we import matrix inverse operator
2. The next step is to define the functions f1, f2, f3 and the Jacobian matrix. For calculating Jacobian matrix, numerical discretization method is used by considering a small grid size.
3. The next step is to define initial conditions for guess values using a column matrix. At the zeroth time step, the new values are equal to the guess values.
4. We initialize the function column matrix by copying from guess matrix to create a new matrix and avoid any errors due to pointer linkage
5. The tolerance and error values are defined along with the step size and relaxation factors
6. The next step is to define a time loop ranging from to zero to length of time array and using multivariate NR method, the new values for Y column matrix is calculated until error becomes less than the tolerance value
7. The values of y1, y2 and y3 are printed for each time step iteration.

The program code is as follows:

''' 
Program to solve the system of equations

dy1/dt = -0.04y1 + 10^4 y2.y3
dy2/dt = 0.04y1 - 10^4 y2.y3 - 3.10^7 y2^2
dy3/dt = 3.10^7 y2^2

Written by Vignesh V
'''

import numpy as np 
from numpy.linalg import inv

#function definition
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)

#function to calculate the Jacobian matrix
def Jacobian(y1_old,y2_old,y3_old,y1,y2,y3,dt):

	J = np.zeros((3,3))
	dy = 1e-4
	
	#Row1
	J[0,0] = (f1(y1_old,y1+dy,y2,y3,dt)-f1(y1_old,y1,y2,y3,dt))/dy	
	J[0,1] = (f1(y1_old,y1,y2+dy,y3,dt)-f1(y1_old,y1,y2,y3,dt))/dy
	J[0,2] = (f1(y1_old,y1,y2,y3+dy,dt)-f1(y1_old,y1,y2,y3,dt))/dy

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

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

	return J

#Initial guess values column matrix
y_old =np.zeros((3,1))

y_old[0]=1
y_old[1]=0
y_old[2]=0

y_new=y_old

#function column matrix
F= np.copy(y_old)


tol= 5e-6
dt = 1e-2
t=np.arange(0,36000,dt)
error = 9e9
alpha = 1
iter = 1

for i in range(0,len(t)):
	while(error>tol):  
		F[0]= f1(y_old[0],y_new[0],y_new[1],y_new[2],dt)
		F[1]= f2(y_old[1],y_new[0],y_new[1],y_new[2],dt)
		F[2]= f3(y_old[2],y_new[0],y_new[1],y_new[2],dt)

		J= Jacobian(y_old[0],y_old[1],y_old[2],y_new[0],y_new[1],y_new[2],dt)

		y_new = y_old - alpha*np.matmul(inv(J),F)
		error = np.max(np.abs(y_new-y_old))
		print('Error = {0}'.format(error))
		y_old = y_new 
	
	print('The values of y1={0} y2={1} y3={2} at iteration ={3}'.format(y_new[0],y_new[1],y_new[2],iter))	
	iter = iter+1
	error=9e9
	

5. Results

The results for time step size of 1e-2 seconds and a tolerance of 5e-6 is

y1 = 0.616,y2=3.01e-06 ,y3=0.821

6. Effect of time step on solution

For an implicit method of iterative solution, the numerical solution is stable at even larger time step sizes. However, the accuracy of the results is poor when calculated with larger step size. This is preferred when processing power and time are major concerns as larger time step can be evaluated in relatively less computational time

The results for varying time step sizes are tabulated below:

Time step size (in secs)	y1	y2	y3	Error	Processing time (in secs)
1.00E-02	0.616	3.01E-06	0.821	2.71E-06	82
1.00E-01	0.2386	1.15E-06	0.825	3.99E-06	18.5
1	0.0904	3.97E-07	0.909	4.67E-06	7.5
10	0.0289	1.22E-07	0.949	4.90E-06	4.3

It can be concluded from the results that with an increase in time step size, the results obtained differ significantly. Thus, to get more accurate results, a smaller time step size should be used.

7. Conclusion

Thus, multivariate Newton Rhapson method is used to solve the coupled ordinary differential equations and the effect of time step size on the accuracy of results is established.