Solving a system of ODEs using Implicit Euler method

In this project, I have discussed and proposed a method to solve a system of stuff ODEs using the first order Implicit Euler method

The ODEs I will be solving are in the form of,

As it can be observed it is a system of coupled nonlinear ODEs,

The solution of this system will explode if we use explicit methods, Hence an implicit formulation has been used

I will be coding the numerical solver for these set of equations using python, and the general logic I have followed will be of this form,

Here we need to use 2 basic loops, 1 for integrating every time step and another to the zero of the nonlinear root using the Newton Rhapson Method

Newton Rhapson formulation

Here for a system of ODEs, the Newton Rhapson equation will be in the following form,

Since we are dealing with a coupled ODE system, we need to use a Jacobian instead of a normal derivative, as the ODEs changes w.r.t all the 3 variables

Now, since we have formulated the Newton Rhapson solver, we need to discuss the input for the solver,

The input for the solver should be the in the form of implicit Euler equation,

This equation needs to be rearranged in the following form to make it solvable, using the Newton Rhapson Algorithm,

Once we rearrange all the 3 equations we will get them in this form,

These will be the input functions for our program,

Here I have pasted the code for my program,

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

"""
State Space Formulation

"""

def S1(y1,y2,y3,y1_new,y2_new,y3_new,h):
	
	return y1 + h*((-0.04*y1_new) + (10000*y2_new*y3_new)) - y1_new
	#return y1 + h*(y2_new) - y1_new

def S2(y1,y2,y3,y1_new,y2_new,y3_new,h):

	return y2 + h*((0.04*y1_new) - (10000*y2_new*y3_new) - (3*10000000*(y2_new**2))) - y2_new
	#return y2 + h*(-0.01*y1_new - 1000*y2_new) - y2_new

def S3(y1,y2,y3,y1_new,y2_new,y3_new,h):

	return y3 + h*(3*10000000*(y2_new**2)) - y3_new
	#return y3 + h*(6*(y2_new)) - y3_new


def jacobian(y1,y2,y3,y1_new,y2_new,y3_new,h):

	J = np.ones((3,3))
	d = 1e-9

	J[0,0] = (S1(y1,y2,y3,(y1_new+d),y2_new,y3_new,h) - S1(y1,y2,y3,y1_new,y2_new,y3_new,h))/d
	J[0,1] = (S1(y1,y2,y3,y1_new,(y2_new+d),y3_new,h) - S1(y1,y2,y3,y1_new,y2_new,y3_new,h))/d
	J[0,2] = (S1(y1,y2,y3,y1_new,y2_new,(y3_new+d),h) - S1(y1,y2,y3,y1_new,y2_new,y3_new,h))/d


	J[1,0] = (S2(y1,y2,y3,(y1_new+d),y2_new,y3_new,h) - S2(y1,y2,y3,y1_new,y2_new,y3_new,h))/d
	J[1,1] = (S2(y1,y2,y3,y1_new,(y2_new+d),y3_new,h) - S2(y1,y2,y3,y1_new,y2_new,y3_new,h))/d
	J[1,2] = (S2(y1,y2,y3,y1_new,y2_new,(y3_new+d),h) - S2(y1,y2,y3,y1_new,y2_new,y3_new,h))/d

	J[2,0] = (S3(y1,y2,y3,(y1_new+d),y2_new,y3_new,h) - S3(y1,y2,y3,y1_new,y2_new,y3_new,h))/d
	J[2,1] = (S3(y1,y2,y3,y1_new,(y2_new+d),y3_new,h) - S3(y1,y2,y3,y1_new,y2_new,y3_new,h))/d
	J[2,2] = (S3(y1,y2,y3,y1_new,y2_new,(y3_new+d),h) - S3(y1,y2,y3,y1_new,y2_new,y3_new,h))/d

	return J

def NewtRhap(y1,y2,y3,y1_guess,y2_guess,y3_guess,h):

	S_old = np.ones((3,1))
	S_old[0] = y1_guess
	S_old[1] = y2_guess
	S_old[2] = y3_guess

	F = np.ones((3,1))
	
	error = 9e9
	tol = 1e-9
	aplha = 1

	iter = 1

	while error > tol:

		J = jacobian(y1,y2,y3,S_old[0],S_old[1],S_old[2],h)

		F[0] = S1(y1,y2,y3,S_old[0],S_old[1],S_old[2],h)
		F[1] = S2(y1,y2,y3,S_old[0],S_old[1],S_old[2],h)
		F[2] = S3(y1,y2,y3,S_old[0],S_old[1],S_old[2],h)

		S_new = S_old - aplha*(np.matmul(inv(J),F))

		error = np.max(np.abs(S_new - S_old))

		S_old = S_new

		iter = iter + 1

	return [S_new[0],S_new[1],S_new[2]] 

def implicit_euler(inty1,inty2,inty3,tspan,dt):
	 
	t = np.arange(0,tspan,dt)
	y1 = np.zeros(len(t))
	y2 = np.zeros(len(t))
	y3 = np.zeros(len(t))

	y1[0] = inty1
	y2[0] = inty2
	y3[0] = inty3

	y1_guess = 10
	y2_guess = 10
	y3_guess = 10
	
	for i in range(1,len(t)):
		
		y1[i] , y2[i] , y3[i] = NewtRhap(y1[i-1],y2[i-1],y3[i-1],y1_guess,y2_guess,y3_guess,dt)
		
		y1_guess = y1[i]
		y2_guess = y2[i]
		y3_guess = y3[i]

	return [t,y1,y2,y3]

t,y1,y2,y3 = implicit_euler(1,0,0,600,0.6)

plt.plot(t,y1,'b')
plt.plot(t,y2,'g')
plt.plot(t,y3,'y')
plt.legend(["graph of y1","graph of y2","graph of y3"])
plt.show()

We call the program with the following initial conditions for y1, y2 and y3 respectively - 1,0,0

We integrate this system for t = 600 seconds and out step size is 0.6

We then get the following results,

From here we can see that graphs y1 and y3 exhibit very stiff behaviour whereas graph 2 hardly changes,

If we had used an explicit solver, our entire solution set would have exploded.