Multivariate Newton Raphson Solver using Python

AIM: To solve the given Stiff ordinary differential equations employing implicit/BDF and solving by Multivariate Newton-Raphson method

`(dy1)/dt = -0.04*y1 + 10^4*y2*y3`

`(dy2)/dt = 0.04*y1 - 10^4*y2*y3 - 3*10^7*y2^2`

`(dy3)/dt = 3*10^7*y2^2`

Employing BDF (backward differencing)/ Implicit method to the ODE\'s we get:

`(y1^n-y1^(n-1))/(deltat) = -0.04*y1 + 10^4*y2*y3`

`(y2^n-y2^(n-1))/(deltat) = 0.04*y1 - 10^4*y2*y3 - 3*10^7*y2^2`

`(y3^n-y3^(n-1))/(deltat) =3*10^7*y2^2`

Rearranging the equations in the form LHS-RHS = 0, we get a set of coupled Non-linear equations:

`f1^n = y1^n + deltat*(-0.04*y1 + 10^4*y2*y3) - y1^(n-1)`

`f2^n = y2^n + deltat*(0.04*y1 - 10^4*y2*y3 - 3*10^7*y2^2) - y2^(n-1)`

`f3^n = y3^n + deltat*(3*10^7*y2^2) - y3^(n-1)`

where

`y^n = y_(New)`

`y^(n-1) = y_(previous)`

y1_guess = 1

y2_guess = 1

y3_guess = 1

`Y_(New) = Y_(old) - alpha*J^-1 * F`

alpha: relaxation factor

F: Function matrix = `[[f_1],[f_2],[f_3]]`

J: Jacobian matrix

`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)]]`

Procedure:

Defining functions and jacobian

Initial guessing of y_previous (to be used in time loop) and y_old (to be used in while loop)

Initiating time loop that would run for n intervals

Initiating while loop

Calculating values of y_new using initial guess values

Updating the old values to again run inside while loop until convergence is achieved

Updating y_previous with y_new obtained after convergence from while loop (to be used for next) time step

Storing y_new values after each time step in a new matrix

Plotting the final values as a function of timesteps after calculations for every time steps are done

Python code:

\"\"\"


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

def f1(y1_previous, y1, y2, y3, dt):
	return y1_previous - 0.04*y1*dt + 10000*y2*y3*dt - y1

def f2(y2_previous, y1, y2, y3, dt):
	return y2_previous + 0.04*y1*dt - 10000*y2*y3*dt - 3*pow(10,7)*pow(y2,2)*dt - y2

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

def jacobian(y1_previous,y2_previous,y3_previous, y1, y2, y3,dt):
	 J = np.ones((3,3))
	 h = 1e-5 #step size

	 # Row 1
	 J[0,0] = (f1(y1_previous, y1+h, y2, y3, dt) - f1(y1_previous, y1, y2, y3, dt))/h
	 J[0,1] = (f1(y1_previous, y1, y2+h, y3, dt) - f1(y1_previous, y1, y2, y3, dt))/h
	 J[0,2] = (f1(y1_previous, y1, y2, y3+h, dt) - f1(y1_previous, y1, y2, y3, dt))/h               

	 # Row 2
	 J[1,0] = (f2(y2_previous, y1+h, y2, y3, dt) - f2(y2_previous, y1, y2, y3, dt))/h
	 J[1,1] = (f2(y2_previous, y1, y2+h, y3, dt) - f2(y2_previous, y1, y2, y3, dt))/h
	 J[1,2] = (f2(y2_previous, y1, y2, y3+h, dt) - f2(y2_previous, y1, y2, y3, dt))/h

	 # Row 3
	 J[2,0] = (f3(y3_previous, y1+h, y2, y3, dt) - f3(y3_previous, y1, y2, y3, dt))/h
	 J[2,1] = (f3(y3_previous, y1, y2+h, y3, dt) - f3(y3_previous, y1, y2, y3, dt))/h
	 J[2,2] = (f3(y3_previous, y1, y2, y3+h, dt) - f3(y3_previous, y1, y2, y3, dt))/h

	 return J

y1_previous = 1
y2_previous = 0
y3_previous = 0

# Defining initial guess
Y_old = np.ones((3,1))
Y_old[0] = y1_previous
Y_old[1] = y2_previous
Y_old[2] = y3_previous

F = np.copy(Y_old)
y_1 = []
y_2 = []
y_3 = []

t = 600
dt = 0.01
n = int(t/dt)
alpha = 1
tol = 1e-6
num_of_time_steps=[]
iter = 1
for i in range (0,n):

	# starting newton raphson loop
	error = 9e9
	
	while error>tol:
		J = jacobian(y1_previous,y2_previous,y3_previous,Y_old[0], Y_old[1], Y_old[2],dt)

		F[0] = f1(y1_previous, Y_old[0], Y_old[1], Y_old[2],dt)
		F[1] = f2(y2_previous,Y_old[0], Y_old[1], Y_old[2],dt)
		F[2] = f3(y3_previous,Y_old[0], Y_old[1], Y_old[2],dt)

		Y_new = Y_old - alpha*np.matmul(inv(J),F)

		error = np.max(np.abs(Y_new - Y_old))

		Y_old = Y_new
	
	log_message = \'iteration # = {0} y1 = {1} y2 = {2} y3 = {3}\'.format(iter, Y_new[0], Y_new[1], Y_new[2])
	print(log_message)
	iter = iter + 1
	time = iter*dt

	y_1.append(Y_new[0])
	y_2.append(Y_new[1])
	y_3.append(Y_new[2])

	#dt= dt + 0.1
	num_of_time_steps.append(time)
	Y_old = Y_new
	y1_previous = Y_new[0]
	y2_previous = Y_new[1]
	y3_previous = Y_new[2] 

plt.plot(num_of_time_steps,y_1)
plt.plot(num_of_time_steps,y_2)
plt.plot(num_of_time_steps,y_3)
plt.xlabel(\'Simulation time\')
plt.ylabel(\'Variable Value\')
plt.grid(\'on\')
plt.legend([\'y1\',\'y2\',\'y3\'])
plt.title(\'Multivariate Newton Raphson method for time step size of 0.01s\')
plt.show()

Results:

y1,y2,y3 as a function of the number of time steps: (time step size: 1s; n=600)

y1,y2,y3 as a function of the number of time steps: (time step size: 0.1s, n=6000)

y1,y2,y3 as a function of the number of time steps: (time step size: 0.01s, n=60000)

Conclusions:

From the results and graphs obtained, it can be seen that there is considerable amount of change in the values of y1 and y3 and it will change further if we simulate it for more than 10 minutes, while the value of y2 remains almost constant for almost complete simualtion. This show the stiffness characteristic off a system where some variables shows a considerable amount of change with respect to time while others show no change

The final values of variable for a time step size of 0.1 s and 0.01 s remain almsot equal, this is due to implicit formualtion of the ODE\'s. Hence to reduce computational time and cost, a time step size of 0.1 s can be finalized for this set coupled equations