SOLVING STIFF ODE SYSTEMS BY USING IMPLICIT EULER METHOD AND A MULTIVARIATE NEWTON RAPHSON SOLVER

           SOLVING STIFF ODE SYSTEMS BY USING IMPLICIT EULER METHOD AND A MULTIVARIATE NEWTON RAPHSON SOLVER

l. OBJECTIVES

1. Write a code in python to solve a system of stiff ODEs using the Implicit Euler Method(Backward Difference Scheme) and the Multivariate Newton Raphson solver.

2. Test the stability and accuracy of the solution at different time steps.

ll. INTRODUCTION

In any combustion mechanism, the rates of production or dissipation of different species can be represented in the form of coupled differential equations. These equations cannot be solved analytically due to the complexity of the equations and hence different solvers are employed to solve them approximately.

1. Newton Raphson Method

In this method, the roots of the function can be determined as follows

• `x^(n+1)=x^n-alpha(f(x^n))/(f'(x^n))`

where the superscript 'n' represents the current timestep.

This method is only applicable when the solution of a single equation needs to be calculated. For a Multivariate system, the roots can be determined as follows-

• `X^(n+1)=X^n-alpha(J^-1F)`

where `Jrarr`Jacobian Matrix,

          `Frarr`Function Matrix

          `alphararr`Relaxation Factor

2. Jacobian Matrix

The Jacobian Matrix represents a matrix of partial derivatives of the different functions with respect to the variables involved in the system.

The number of rows and the column of the matrix is equal to the number of functions and variables involved in the given system.

For a system having three functions and three variables, the Jacobian matrix can be represented as follows-

• `J=[[(delf_1)/(delx_1),(delf_1)/(delx_2),(delf_1)/(delx_3)],[(delf_2)/(delx_1),(delf_2)/(delx_2),(delf_2)/(delx_3)],[(delf_3)/(delx_1),(delf_3)/(delx_2),(delf_3)/(delx_3)]]`

lll. PROBLEM STATEMENT

1. Solve the following system of coupled ODEs-

• `(dy_1)/(dt)=-0.04y_1+10^4y_2y_3`
• `(dy_2)/(dt)=0.04y_1-10^4y_2y_3-3.10^7y_2^2`
• `(dy_3)/(dt) = 3.10^7y_2^2`

2. Solve using the Implicit Euler Method (Backward Differencing Scheme) and the multivariate Newton Raphson Solver.

3. Simulate for 10 minutes.

lV. ASSUMPTIONS

1. The initial guess values are the three variables are-

• `y_1=1,y_2=0,y_3=0`

2. The relaxation factor is assumed as 1

3. The tolerance of the solution is assumed to be 1e-10.

V. METHODOLOGY

1. Defining the Function Matrix

To solve the given equations, the backward differencing scheme is employed -

• `(y_1^n-y_1^(n-1))/(Deltat)=-0.04y_1^n+10^4y_2^ny_3^n`
• `(y_2^n-y_2^(n-1))/(Deltat)= 0.04y_1^n-10^4y_2^ny_3^n-3.10^7(y_2^n)^2`
• `(y_3^n-y_3^(n-1))/(Deltat)=3.10^7(y_2^n)^2`    

substituting `h=Deltat`and rearranging the given equations -

• `y_1^n=y_1^(n-1)+h(-0.04y_1^n+10^4y_2^ny_3^n)`
• `y_2^n=y_2^(n-1)+h(0.04y_1^n-10^4y_2^ny_3^n-3.10^7(y_2^n)^2)`
• `y_3^n=y_3^(n-1)+h(3.10^7(y_2^n)^2)`

using the above equations, the function vector F can be represented as-

• `f_1=y_1^(n-1)+h(-0.04y_1^n+10^4y_2^ny_3^n)-y_1^n`
• `f_2=y_2^(n-1)+h(0.04y_1^n-10^4y_2^ny_3^n-3.10^7(y_2^n)^2)-y_2^n`
• `f_3=3.10^7(y_2^n)^2-y_3^n`

The Function matrix is given as follows-

• `F=[[f_1],[f_2],[f_3]]`

2. Define the Jacobian Matrix

As the jacobian matrix represents a matrix of partial derivatives of the different functions with respect to the variables involved in the system, we need to calculate the partial derivative of the given equations.

To calculate the partial derivatives, we shall use the following approximations -

• `(delf)/(dely)=(f(y+dy)-f(y))/(dy)`,where `dely=1e-6`

The jacobian natrix for the given system of eqautions is represented as follws-

• `J=[[(delf_1)/(dely_1),(delf_1)/(dely_2),(delf_1)/(dely_3)],[(delf_2)/(dely_1),(delf_2)/(dely_2),(delf_2)/(dely_3)],[(delf_3)/(dely_1),(delf_3)/(dely_2),(delf_3)/(dely_3)]]`

3. Define the Initial guess values

The initial guess values of the three variables are -

• `y_1=1,y_2=0,y_3=0`

4. Define the Time and Timestep

The start and end time is as given in the problem statement-

• t_start = 0 seconds and t_end = 10 minutes = 600 seconds

Time step size h= 0.6

5. Initialize a solution array

A solution array is initialized to store the solution at the end of every time step.

6. Define the Outer Loop

A column vector is initialized each for the value of the variables at the previous time step, the guess value of the variable, and the function vectors.

The parameters for the Newton Raphson solver such as the relaxation factor, tolerance, and error are set and the inner loop is defined.

7. Define the Inner Loop

The jacobian and Function vectors are computed and then the solution is computed using the multivariate Newton Raphson method as follows-

•  `Y^(n+1)=Y^n-alpha(J^-1F)`

The error is computed and the guess values of the variables are updated. The loop repeats until the error becomes less than the tolerance. Finally, the results are stored and the time step is updated. The process continues and the values of the variables are computed at each time step.

8. Plot the Graph

The graph of the solutions of the variables is plotted against  the time steps

Vl. PYTHON CODE

"""

	Program to solve the system of Stiff ODEs

	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

"""

from numpy import*
from numpy.linalg import inv
from matplotlib.pyplot import*

a = 0.04
b = 1e4
c = 3e7
# Defining the Non-Linear Equation Function
def f1(y10,y1, y2, y3, h):
	return (y10 + h*(-a*y1 + b*y2*y3) - y1)

def f2(y20,y1, y2, y3, h):	 
	return (y20 + h*(a*y1 - (b*y2*y3) - (c*y2*y2)) - y2)

def f3(y30,y1, y2, y3, h):
	return (y30 + h*(c*y2*y2) - y3)	

# Defining the Jacobian Function	

def jacobian(y10,y20,y30,y1,y2,y3,h):

	J = ones((3,3))
	dy = 1e-6
	y1_prime = y1 + dy
	y2_prime = y2 + dy
	y3_prime = y3 + dy

	# Row 1
	J[0,0] = (f1(y10, y1_prime, y2, y3, h) - f1(y10, y1, y2, y3, h)) / dy
	J[0,1] = (f1(y10, y1, y2_prime, y3, h) - f1(y10, y1, y2, y3, h)) / dy
	J[0,2] = (f1(y10, y1, y2, y3_prime, h) - f1(y10, y1, y2, y3, h)) / dy

	# Row 2
	J[1,0] = (f2(y20, y1_prime, y2, y3, h) - f2(y20, y1, y2, y3, h)) / dy
	J[1,1] = (f2(y20, y1, y2_prime, y3, h) - f2(y20, y1, y2, y3, h)) / dy
	J[1,2] = (f2(y20, y1, y2, y3_prime, h) - f2(y20, y1, y2, y3, h)) / dy

	# Row 3
	J[2,0] = (f3(y30, y1_prime, y2, y3, h) - f3(y30, y1, y2, y3, h)) / dy
	J[2,1] = (f3(y30, y1, y2_prime, y3, h) - f3(y30, y1, y2, y3, h)) / dy
	J[2,2] = (f3(y30, y1, y2, y3_prime, h) - f3(y30, y1, y2, y3, h)) / dy

	return J
# Defining the initial guess values
y10 = 1
y20 = 0
y30 = 0
# Column Vector
y_old = ones((3,1))
y_old[0] = y10
y_old[1] = y20
y_old[2] = y30
# Column vector for F
F = copy(y_old)
# Defining the Time and Timestep	
h = 0.1
t = 600
n=int(t/h)
# Parameter for Newton Raphson Solver
alpha = 1
tol = 9e-6
iter = 1
# Initial error
error = 9e9
# Final Values
x1 = []
x2 = []
x3 = []
iteration = []
#Starting the newton rapson loop
iter = 1
# Outer Time Loop
for i in range (0,n):
	
	while error > tol:
	    # Computing the Jacobian
		J = jacobian(y10,y20,y30,y_old[0],y_old[1],y_old[2], h)
	    # Computing the F Vector
		F[0] = f1(y10, y_old[0], y_old[1], y_old[2], h)
		F[1] = f2(y20, y_old[0], y_old[1], y_old[2], h)
		F[2] = f3(y30, y_old[0], y_old[1], y_old[2], h)
	    # Computing New Values
		y_new = y_old - alpha*matmul(inv(J), F)
	    # Computing the maximum absolute error
		error = max(abs(y_new-y_old))
		print(error)
	    # Updating old values
		y_old =y_new
	#log message	
	log_message = 'iteration No. ={0} y1 = {1} y2 = {2} y3 = {3}'.format(iter,y_new[0],y_new[1],y_new[2])
	print(log_message)

	#storing the value for plotting
	x1.append(y_new[0])
	x2.append(y_new[1])
	x3.append(y_new[2])
	iter = iter+1
	iteration.append(iter)
	y_old = y_new
	y10 = y_new[0]
	y20 = y_new[1]
	y30 = y_new[2]
	error = 9e9

#Plotting the results

plot(iteration,x1,color = 'b')
plot(iteration,x2,color = 'r')
plot(iteration,x3,color = 'k')
xlabel('Number of steps')
ylabel('Variables(species)')
legend(['y1','y2','y3'])
title('Solution of ODE')
show()	

Vll. OUTPUTS

Vlll. RESULTS

1. It can be observed from the plots that as the time step size decreases, the accuracy of the solution increases. However, the time required to compute the solution also increases.

2. The solution is stable even when the solver is given a higher time step. This shows that the implicit solver is highly stable and a solution can be obtained, although the solution may not be accurate.

3. The determination of the size of the time step will be a compromise between accuracy and cost depending on the field and scope of application.

4. The ideal 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 time step which would give the required accuracy will be the ideal time step.

lX. CONCLUSION

A system of stiff ODEs is solved using the Implicit Euler Method and the Multivariate Newton Raphson solver and the stability and accuracy of the solution are tested at different time steps.