Week 6 - Multivariate Newton Rhapson Solver

Calculating the solution to the non linear coupled ODEs by multivariate newton raphson method

Explanation of the method and approach

The derivative in the ODEs can be written by backward differencing i.e.

(y_1(n)-y_1(n-1))/dt=-0.04*y_1+(10)^4*y_2*y_3

(y_2(n)-y_2(n-1))/dt=0.04*y_1-(10)^4*y_2*y_3-3*(10)^7*(y_2)^2

(y_3(n)-y_3(n-1))/dt=3*(10)^7*(y_2)^2

If every thing is taken to the left hand side of the equation, then the ODEs can be written as:

Equation 1-----(y_1(n)-y_1(n-1))-dt*(-0.04*y_1+(10)^4*y_2*y_3)=0

Equation2------(y_2(n)-y_2(n-1))-dt*(0.04*y_1-(10)^4*y_2*y_3-3*(10)^7*(y_2)^2)=0

Equation3------(y_3(n)-y_3(n-1))-dt*(3*(10)^7*(y_2)^2)=0

Here, y(n-1) is the value at old time step, initially we assumed a value (given in the question, the value y at time step 0)of

y_0=1

y_1=0

y_2=0

This is the value at time step 0, running the iterations for 10 min i.e 600 sec

Taking dt =0.01

total 600/0.01=60000 no.of iterations will be run.

The multivariate newton raphson method is:

y(n)=y(n-1)-alpha*(inv(J)*F)

alpha=1

J=jacobian matrix=

df_1/dy_1 df_1/dy_2 df_1/dy_3

df_2/dy_1 df_2/dy_2 df_2/dy_3

df_3/dy_1 df_3/dy_2 df_3/dy_3

The jacobian matrix calculated numerically will give,

h=step increment taken for y_1, y_2 , y_3 =1e-4

doing forward differencing here,

df_1/dy_1=(f_1(y_1+h,y_2,y_3)-f_1(y_1,y_2,y_3))/h

df_1/dy_2=(f_1(y_1,y_2+h,y_3)-f_1(y_1,y_2,y_3))/h

df_1/dy_3=(f_1(y_1,y_2,y_3+h)-f_1(y_1,y_2,y_3))/h

df_2/dy_1=(f_2(y_1+h,y_2,y_3)-f_2(y_1,y_2,y_3))/h

df_2/dy_2=(f_2(y_1,y_2+h,y_3)-f_2(y_1,y_2,y_3))/h

df_2/dy_3=(f_2(y_1,y_2,y_3+h)-f_2(y_1,y_2,y_3))/h

df_3/dy_1=(f_3(y_1+h,y_2,y_3)-f_3(y_1,y_2,y_3))/h

df_3/dy_2=(f_3(y_1,y_2+h,y_3)-f_3(y_1,y_2,y_3))/h

df_3/dy_3=(f_3(y_1,y_2,y_3+h)-f_3(y_1,y_2,y_3))/h

We also assume guess values of y_1, y_2,y_3 to start performing the newton raphson technique iterations, here they are taken as

y1_guess=1

y2_guess=0

y3_guess=0

F matrix is a 3 x 1 matrix and

F[0]= value of function f_1 with guess values

F[1]= value of function f_2 with guess values

F[2]=value of function f_3 with guess values

The value at time step 1 is calculated by the assuming some initial guess values,

y(value at new time step)=y(value assumed)-alpha*inv(J(for assumed value))*F(at assumed value)

This iteration will be run till the error is greater than tolerance, and everytime the y_guess values are replaced with the new values of y.

This same process will be done for each time step, the final new value of y after each iteration that we get will be taken as value at that time step and for the next time step,it will act as the old time step value.

We guess a value everytime in the before starting the newton raphson iteration in every time step.

error=9e9

tolerance =1e-12

The code:

#final time loop
import numpy as np
from numpy.linalg import inv
'''
odes to solve are

4*x_1*x_1-625*x_2*x_2+2*x_3-1=0
20*x_3+np.exp(-x_1*x_2)+9=0


y_1=1
y_2=0
y_3=0
'''
Y_old=np.ones((3,1))
Y_old[0]=1
Y_old[1]=0
Y_old[2]=0
dt=0.01
h=1e-4
for t in np.arange(0,600,dt):
	
	def f_1(y_1,y_2,y_3):
		return y_1-Y_old[0]-dt*(-0.04*y_1+pow(10,4)*y_2*y_3)
	def f_2(y_1,y_2,y_3):
		return y_2-Y_old[1]-dt*(0.04*y_1-pow(10,4)*y_2*y_3-3*pow(10,7)*pow(y_2,2))
	def f_3(y_1,y_2,y_3):
		return y_3-Y_old[2]-dt*(3*pow(10,7)*pow(y_2,2))
	#defining jacobian
	def jacobian(y_1,y_2,y_3):
		J=  np.ones((3,3))
		#derivative of f_1
		J[0,0]=(f_1(y_1+h,y_2,y_3)-f_1(y_1,y_2,y_3))/h
		J[0,1]=(f_1(y_1,y_2+h,y_3)-f_1(y_1,y_2,y_3))/h
		J[0,2]=(f_1(y_1,y_2,y_3+h)-f_1(y_1,y_2,y_3))/h
		#derivative of f_2
		J[1,0]=(f_2(y_1+h,y_2,y_3)-f_2(y_1,y_2,y_3))/h
		J[1,1]=(f_2(y_1,y_2+h,y_3)-f_2(y_1,y_2,y_3))/h
		J[1,2]=(f_2(y_1,y_2,y_3+h)-f_2(y_1,y_2,y_3))/h
		#derivative of f_3
		J[2,0]=(f_3(y_1+h,y_2,y_3)-f_3(y_1,y_2,y_3))/h
		J[2,1]=(f_3(y_1,y_2+h,y_3)-f_3(y_1,y_2,y_3))/h
		J[2,2]=(f_3(y_1,y_2,y_3+h)-f_3(y_1,y_2,y_3))/h
		return J


	# we copied the 3 x 1 matrix of x_old to just save some lines of creating new np.ones((3,1)) column matrix for F matrix, F is totally different
	#from x_old matrix, we could have also defined a different matrix as F=np.ones((3,1)) for F, because anyhow we are changing its value and re
	#defining the value of elements of F inside the loop, we just need a skeleton 3 x 1 matrix having 1 as each element to have a matrix'body'
	
	error=9e9
	
	F=np.ones((3,1))
	tol=1e-12
	alpha=1
	iter=0
	Y_guess=np.ones((3,1))
	Y_guess[0]=1
	Y_guess[1]=0
	Y_guess[2]=0
	while error > tol:
		J=jacobian(Y_guess[0],Y_guess[1],Y_guess[2]);
		F[0]=f_1(Y_guess[0],Y_guess[1],Y_guess[2]);
		F[1]=f_2(Y_guess[0],Y_guess[1],Y_guess[2]);
		F[2]=f_3(Y_guess[0],Y_guess[1],Y_guess[2]);
		Y_new=Y_guess-alpha*np.matmul(inv(J), F)
		error=np.max(np.abs(Y_new-Y_guess))
		
		Y_guess=Y_new
		iter=iter+1
		#log_message='iteration # = {0} y_1 ={1} y_2 = {2} y_3 ={3}'.format(iter,Y_new[0],Y_new[1],Y_new[2])
		#print(log_message)
	#print(error)	
	#print(iter)
	Y_old=Y_new
	print(t)
	log_message='iteration # = {0} y_1 ={1} y_2 = {2} y_3 ={3}'.format(iter,Y_new[0],Y_new[1],Y_new[2])
	print(log_message)
	#print(f_1(Y_new[0],Y_new[1],Y_new[2]))
	y1_final=Y_new[0]
	y2_final=Y_new[1]
	y3_final=Y_new[2]
	
	#print(y1_final)
	#print(y2_final)
	#print(y3_final)
	#print(error)
print('finish')

Observation:

The last five iterations are given here,

599.95
iteration # = 19 y_1 =[0.39993503] y_2 = [2.63134078e-06] y_3 =[0.60006234]
5.260236690673992e-13
599.96
iteration # = 19 y_1 =[0.39993296] y_2 = [2.63131842e-06] y_3 =[0.60006441]
5.259126467649367e-13
599.97
iteration # = 19 y_1 =[0.39993088] y_2 = [2.63129605e-06] y_3 =[0.60006649]
5.259126467649367e-13
599.98
iteration # = 19 y_1 =[0.3999288] y_2 = [2.63127369e-06] y_3 =[0.60006857]
5.259126467649367e-13
599.99
iteration # = 19 y_1 =[0.39992672] y_2 = [2.63125132e-06] y_3 =[0.60007064]

The last iteration will be, after end of 10 minutes.

y_1=0.39992672

y_2=2.63125132e-06

y_3=0.60007064

As this is a implicit solution, the dt value does not have a restriction because the answer is always stable but keeping the accuracy in mind , we should take low values of time step . Here we take time step as 0.01 , but taking the time step as 0.1 also gives us the same answer, i.e taking a bigger time step also gives the accurate answer here, so the taken time step works as ideal time step here.

The answer comes at time step 0.1 comes as,

y_1=0.39997575

y_2=2.63177929e-06

y_3=0.6000216

which is the same as that when time step is 0.01, so the accuracy is maintained even at 0.1 time step.