Multivariate Newton Raphson

Problem

• Write a Python program to solve the problem using Implicit Euler Method/Backward Differencing.
• Incorporate numerical jacobian which should be estimated numerically and not analytically.
• Run the simulation for 10 minutes.
• Choose a suitable step size and explain the reason for choosing that step size.

ODE system

`(dyunderset(1)())/dt =-0.04yunderset(1)() +10^4yunderset(2)()yunderset(3)()`

`(dyunderset(2)())/dt=0.04yunderset(1)()-10^4yunderset(2)()yunderset(3)()-3*10^7yunderset(2)()^2`

`(dyunderset(3)())/dt=3*10^7yunderset(2)()^2`

Newton Raphsom Method

Newton Raphson Method is an iterative method which is used to find roots. We used the implicit differential method because it is a stiff system.

`Yunderset(n+1,New)()=Yunderset(n,old)()-alpha*J^-1*F`

`Yunderset(n+1,New)()`=new value updated after every iteration

`Yunderset(n,old)()`= old value

`alpha`=relaxation factor

`J^-1`= Jacobian matrix

Jacobian matrix has been calculates numerically

`J=[[(partialfunderset(1)())/(partialyunderset(1)()),(partialfunderset(1)())/(partialyunderset(2)()),(partialfunderset(1)())/(partialyunderset(3)())],[(partialfunderset(2)())/(partialyunderset(1)()),(partialfunderset(2)())/(partialyunderset(2)()),(partialfunderset(2)())/(partialyunderset(3)())],[(partialfunderset(3)())/(partialyunderset(1)()),(partialfunderset(3)())/(partialyunderset(2)()),(partialfunderset(3)())/(partialyunderset(3)())]]`

were

`(df)/(dy)=(f(x-h)-f(x))/h`

and

`F`= Function matrix

`funderset(1)()=yunderset(1)()^(n-1)+h(-0.04yunderset(1)()^(n)+10^4yunderset(2)()^nyunderset(3)()^n)-yunderset(1)()^n=0`

`funderset(2)()=yunderset(2)()^(n-1)+h(0.04yunderset(1)()^(n)-10^4yunderset(2)()^nyunderset(3)()^n-3*10^7yunderset(2)()^2)-yunderset(2)()^n=0`

`funderset(3)()=yunderset(2)()^(n-1)+h(3*10^7yunderset(3)()^2)-yunderset(2)()^n=0`

Initial values

y1=1

y2=0

y3=0

Python code

\"\"\"
Program to solve a differential non-linear coupled equations using th miltivariet Newton Raphson  Method
	ODE System
	dy1/dt=-0.04y1+pow(10,4)*y2y3
	dy2/dt=0.04y1-pow(10,4)*y2y3-3*pow(10,7)*pow(y2,2)
	dy3/dt=3*pow(10,7)*pow(y2,2)

	#BDF using Multivariate Newton Raphson Method

	f1=y1(n-1) + dt*(-0.04y1(n) +pow(10,4)*y2(n)y3(n)) -y1(n) =0
	f2=y2(n-1) + dt*(0.04y1(n) -pow(10,4)*y2(n)y3(n) -3*pow(10,7)*pow(y2(n),2)) -y2(n) =0
	f3=y3(n-1) + dt*(3*pow(10,7)*pow(y2(n),2)) -y3(n) =0
\"\"\"

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


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

def jacobian (y1_old, y2_old, y3_old, y1, y2, y3, dt):

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

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

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

	return J

#Initial values	
y1 = 1
y2 = 0
y3 = 0

#Time loop
dt = 1#step
t = np.arange(0, 120000,dt)

#Final values
Y1_new=[]
Y2_new=[]
Y3_new=[]

iteration=0
for i in range(1,len(t)+1):

	#parameters for the newton rhapson solver
	error = 9e9
	tolerance = 1e-8
	alpha = 1
	
	#Old values
	
	Y_old = np.ones((3,1))
	Y_old[0] = y1
	Y_old[1] = y2
	Y_old[2] = y3

	#guess values	
	Y_guess = np.ones((3,1))
	Y_guess[0] = y1
	Y_guess[1] = y2
	Y_guess[2] = y3

	F=np.copy(Y_guess)

	while error>tolerance:

		J = jacobian(Y_old[0], Y_old[1], Y_old[2], Y_guess[0], Y_guess[1], Y_guess[2], dt)
		
		F[0] = f1(Y_old[0], Y_guess[0], Y_guess[1], Y_guess[2], dt)
		F[1] = f2(Y_old[1], Y_guess[0], Y_guess[1], Y_guess[2], dt)
		F[2] = f3(Y_old[2], Y_guess[0], Y_guess[1], Y_guess[2], dt)
		
		Y_new = Y_guess - alpha*np.matmul(inv(J),F)

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

		Y_guess = Y_new

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

		iteration = iteration + 1

	y1 = Y_new[0]
	y2 = Y_new[1]
	y3 = Y_new[2]

	Y1_new.append(Y_new[0])
	Y2_new.append(Y_new[1])
	Y3_new.append(Y_new[2])

plt.plot(t,Y1_new,color=\'purple\')
plt.plot(t,Y2_new,color=\'blue\')
plt.plot(t,Y3_new,color=\'pink\')
plt.xlabel(\'Time [s]\')
plt.legend([\'Y1\',\'Y2\',\'Y3\'])
plt.show()

Result

Observation

An implicit fuction is not sesitive to time step, we used different time steps lesser than 1 and we could not see a difference between the results. Changing the time step just helps to increase or decrease the real time of the simulation without change the results showed on the graph.
A total period of 10 minutes in real time are equal to inner loop 120000 with a time step equal to 1.