Simulating a simple stiff non-linear coupled ODE system using Newton-Raphson solver with implicit time integration

objective:

To simulate a simple stiff non-linear coupled 3-ODE system using the multi-variate Newton- Raphson solver with implicit time integration for a time of 10 minutes. The ODEs are as follows;

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

`dy_2/dt = 0.004 y_1 - 10^4 y_2y_3 - 3*10^7 y_2^2`

`dy_3/dt = 3*10^7 y_2^2`

Solution:

Since the given ODEs are non-linear, they cannot be solved directly like how we usually solve the linear systems. However, they can be converted into root finding problems and solved using Newton-Raphson solver. And since the ODEs are stiff, they can easily become unstable and require very small time step size, which can become computationally costly as our simulation time is large. Therefore, we use the implicit method for time integration as it is unconditionally stable i.e., for any time step size, though compromising a small bit on accuracy.

In order to solve the given stiff non-linear ODE system, we must first discretize the given ODEs implicitly. Implicit time integration means that we\'re discretizing the time derivative using backward difference.

Discretizing the first equation, we get

`(y_1^(t) - y_1^(t-Delta t))/(Delta t) = -0.004 y_1^t + 10^4 y_2^t y_3^t`

Rearranging terms,

`y_1^(t) = y_1^(t-Delta t) + Delta t(-0.004 y_1^t + 10^4 y_2^t y_3^t)`

Similarly, the other 2 equations after discretization can be written as

`y_2^t = y_2^(t-Delta t) + Delta t(0.004 y_1^t - 10^4 y_2^t y_3^t - 3*10^7 (y_2^t)^2)`

`y_3^t = y_3^(t-Delta t) + Delta t (3*10^7 (y_2^t)^2)`

Converting the above discretized equations into root finding problems, we get

`y_1^(t-Delta t) + Delta t(-0.004 y_1^t + 10^4 y_2^t y_3^t) - y_1^(t) = 0`

`y_2^(t-Delta t) + Delta t(0.004 y_1^t - 10^4 y_2^t y_3^t - 3*10^7 (y_2^t)^2) - y_2^t = 0`

`y_3^(t-Delta t) + Delta t (3*10^7 (y_2^t)^2) - y_3^t = 0`

Newton-Raphson solver

The Newton-Raphson solver is used to find the roots of a non-linear equation or a polynomial. It is an iterative method and the solution can be determined with the desired level of accuracy. A simple Newton-Raphson iterative solution is given by;

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

where,

`y^n` - the old guess value

`y^(n+1)` - the next guess value

`alpha` - the relaxation factor

`f(x)` - root finding problem function

`f\'(x)` - first order derivative of f(x) with respect to y

But, since our system has 3 ODEs and 3 unknowns, the simple form of the Newton-Raphson method is not enough. We have to use the multivariate Newton-Raphson solver.

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

where

`[Y]^n` =

`[Y]^(n+1)` =

`J` = Jacobian matrix =

`F` = function vector =

The Jacobian is a matrix of first order derivatives of all functions with respect to all individual unknowns. The first order derivative can be computed either analytically or numerically. In this project, since we\'re required to compute numerically, I\'m using forward difference approximation.

For a function f, the forward difference approximation for the first-order derivative with respect to a variable y is given by

`(del f)/(del y) = (f(y+Delta y) - f(y))/(Delta y)`

Assuming a small value for dt, we can calculate the derivatives for given values of y.

Initial conditions

Since this is an initial value problem, we need to specify the initial conditions i.e., the values of y1, y2, and y3 at time t=0. The initial condition is given as

y1 = 1

y2 = 0

y3 = 0

Time step size

The simulation is to be run for 10 minutes i.e., 600 seconds, and since we are using an implicit time integration scheme, which is unconditionally stable, the time step can be chosen to be 1 second. A smaller step size means more computational time, and a bigger step size means less accuracy. Therefore, it is required to choose a time step that balances both depending on the level of accuracy required.

Ideal time step size?

An implicit scheme as we know is unconditionally stable which means it is stable for all values of time-step size. Therefore, as far as just stability is concerned, there wouldn\'t be any ideal time step size for an implicit scheme. However, there are two factors that must be weighed while considering the time step size for an implicit scheme, the accuracy and computational time. If the time step is too small it might take a long time to solve the system, and if it is too large the inaccuracy in the solution might become significant, all while keeping the solution stable. Therefore, an ideal time step would be one that can achieve the desired level of accuracy in the least amount of computational time.

Program:

\"\"\"
Solving a stiff coupled non-linear 3-ODE system using implicit time integration

dy1/dt = -0.004 y1 + 1e4 y2y3
dy2/dt = 0.004 y1 - 1e4 y2y3 - 3e7 y2^2
dy3/dt = 3e7 y2^2
Discretizing the derivatives using backward difference (Implicit)

y1_(n) = y1_(n-1) + dt (-0.04 y1_(n) + 1e4 y2_(n) y3_(n))
y2_(n) = y2_(n-1) + dt (0.04 y1_(n) - 1e4 y2_(n) y3_(n) - 3e7 (y2_(n))^2)
y3_(n) = y3_(n-1) + dt (3e7 (y2_(n))^2)

Converting into root finding problem functions

y1_(n-1) + dt (-0.04 y1_(n) + 1e4 y2_(n) y3_(n)) - y1_(n) = 0
y2_(n-1) + dt (0.04 y1_(n) - 1e4 y2_(n) y3_(n) - 3e7 (y2_(n))^2) - y2_(n) = 0
y3_(n-1) + dt (3e7 (y2_(n))^2) - y3_(n) = 0
\"\"\"

import numpy as np
from numpy.linalg import inv
import matplotlib.pyplot as plt

def f1(y1_prev, y1, y2, y3, dt):

	return y1_prev + dt*(-0.04*y1 + 1e4*y2*y3) - y1

def f2(y2_prev, y1, y2, y3, dt):

	return y2_prev + dt*(0.04*y1 - 1e4*y2*y3 - 3e7*y2**2) - y2

def f3(y3_prev, y1, y2, y3, dt):

	return y3_prev + dt*(3e7*y2**2) - y3

def jac(Y_prev, Y, dy, dt):

	J = np.zeros((3,3))

	# 1st row
	J[0,0] = (f1(Y_prev[0], Y[0]+dy, Y[1], Y[2], dt) - f1(Y_prev[0], Y[0], Y[1], Y[2], dt))/dy
	J[0,1] = (f1(Y_prev[0], Y[0], Y[1]+dy, Y[2], dt) - f1(Y_prev[0], Y[0], Y[1], Y[2], dt))/dy
	J[0,2] = (f1(Y_prev[0], Y[0], Y[1], Y[2]+dy, dt) - f1(Y_prev[0], Y[0], Y[1], Y[2], dt))/dy

	# 2nd row
	J[1,0] = (f2(Y_prev[0], Y[0]+dy, Y[1], Y[2], dt) - f2(Y_prev[0], Y[0], Y[1], Y[2], dt))/dy
	J[1,1] = (f2(Y_prev[0], Y[0], Y[1]+dy, Y[2], dt) - f2(Y_prev[0], Y[0], Y[1], Y[2], dt))/dy
	J[1,2] = (f2(Y_prev[0], Y[0], Y[1], Y[2]+dy, dt) - f2(Y_prev[0], Y[0], Y[1], Y[2], dt))/dy

	# 3rd row
	J[2,0] = (f3(Y_prev[0], Y[0]+dy, Y[1], Y[2], dt) - f3(Y_prev[0], Y[0], Y[1], Y[2], dt))/dy
	J[2,1] = (f3(Y_prev[0], Y[0], Y[1]+dy, Y[2], dt) - f3(Y_prev[0], Y[0], Y[1], Y[2], dt))/dy
	J[2,2] = (f3(Y_prev[0], Y[0], Y[1], Y[2]+dy, dt) - f3(Y_prev[0], Y[0], Y[1], Y[2], dt))/dy

	return J

# initial guess values
y1 = 1
y2 = 0
y3 = 0

# creating a copy of initial values
y1_copy = [y1]
y2_copy = [y2]
y3_copy = [y3]

# step control
t_end = 600
dt = 2
t = np.arange(0, t_end+dt, dt)
dy = 1e-4

# initial condition
Y = np.zeros((3,1))
Y[0] = y1
Y[1] = y2
Y[2] = y3
Y_prev = np.copy(Y)
Y_copy = []
Y_copy.append(Y)

# parameters for Newton-Raphson solver
tol = 1e-16
alpha = 1

for i in range(1,len(t)):

	# initializing vectors
	Y_old = np.zeros((3,1))
	F = np.copy(Y_old)

	# initial error
	err = 1

	# Newton-Raphson loop
	while err>tol:

		# computing the function vector
		F[0] = f1(Y_prev[0], Y[0], Y[1], Y[2], dt)
		F[1] = f2(Y_prev[1], Y[0], Y[1], Y[2], dt)
		F[2] = f3(Y_prev[2], Y[0], Y[1], Y[2], dt)

		# computing the Jacobian matrix
		J = jac(Y_prev, Y, dy, dt)

		# calculating the new values
		Y = Y_old - alpha*np.matmul(inv(J),F)

		# calculating error
		err = np.max(np.abs(Y - Y_old))

		# updating old values
		Y_old = np.copy(Y)

	# updating previous values
	Y_prev = np.copy(Y)
	y1_copy.append(Y[0])
	y2_copy.append(Y[1])
	y3_copy.append(Y[2])

print(Y)	
plt.plot(t, y1_copy, label = \'y1\')
plt.plot(t, y2_copy, label = \'y2\')
plt.plot(t, y3_copy, label = \'y3\')
plt.title(\'Simulation of a stiff non-linear coupled 3-ODE system for 10 minutes\\nusing Newton-Raphson method with implicit time integration\')
plt.xlabel(\'Time (s)\')
plt.ylabel(\'Magnitude\')
plt.legend(loc = \'upper right\')
plt.show()

Result:

From the above graph, we can see the variation of the quantities y1, y2, and y3 as governed by the given ODE system from t = 0 through t = 10 minutes. This system mimics a chemical reaction and the above graph can be interpreted as the variation of concentrations of species y1, y2, and y3 with time in the chemical system.