Week 6 - Multivariate Newton Rhapson Solver
#Aim-To solve the set of given Ordinary Differential Equations.
`dy_1/(dt)=-0.04y_1+10^4*y_2*y_3`
`dy_2/(dt)=0.04y_1-10^4y_2y_3-3*10^7y_2^2`
`dy_3/(dt)=3*10^7y_2^2`
f1= (y1)(n)-(y1)(n-1)-dt*(-0.04*y1(n)+10^4*y2(n)*y3(n)) =0
f2= (y2)(n)-(y2)(n-1)-dt*(0.04*(y1)(n)-10^4*(y2)(n)*(y3)(n)-3*10^7*((y2)(n))^2 =0
f3= (y3)(n)-(y3)(n-1)-dt*(3*10^7*((y2)(n))^2) =0
For getting value of `Y_(guess)` at next iteration which is `Y_(n ew)`,
`Y_(n ew)=Y_(guess)-alphaJ^-1F`
`Y_(n ew)=[[y_(1(n ew))],[y_(2(n ew))],[y_(3(n ew))]]`
`Y_(guess)=[[y_(1(guess))],[y_(2(guess))],[y_(3(guess))]]`
alpha=Relaxation parameter.
`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)]]`
`F=[[delf1],[delf2],[delf3]]`
#Introduction:This challenge is in resemblence with chemical reaction equations.Like here we are solving for change in value of y1,y2,y3 with respect to time,in their places there could be different species participating in a reaction.So by solving these equations we can know how a particular species involved in a reaction is changing with time like if combustion is taking place ,so concentration of species will definitely change as reaction proceed with time.
#Coding Approach:
1.Required modules were imported.
2.Three functions f1,f2,f3 defined ,solving given three equations.
3.Jacobian function defined for calculating jacobian matrix.
4.Initial condition given for y1,y2,y3 used.
5.Total time for solution kept 10 min=600sec and time step size =0.01 sec.
6.There is one outer time marching loop and inside that is while loop for setting convergence criteria of solution variables.
7.After getting solution ,plot it with respect to time.
import numpy as np
from numpy.linalg import inv
import math
import matplotlib.pyplot as plt
"""
f1= (y1)(n)-(y1)(n-1)-dt*(-0.04*y1(n)+10^4*y2(n)*y3(n))=0
f2= (y2)(n)-(y2)(n-1)-dt*(0.04*(y1)(n)-10^4*(y2)(n)*(y3)(n)-3*10^7*((y2)(n))^2=0
f3=(y3)(n)-(y3)(n-1)-dt*(3*10^7*((y2)(n))^2)=0
"""
m=0.04 #Constants are defined in terms of m,n,p
n=10**4
p=3*(10**7)
def f1(y1_old,y1,y2,y3,dt):
"""
Function takes values of y1,y1_old,y2,y3 and dt and returning value which is going to become zero as we approach correct solution
"""
return y1-y1_old-dt*(-m*y1+n*y2*y3)
def f2(y2_old,y1,y2,y3,dt):
"""
Function takes values of y1,y2_old,y2,y3 and dt and returning value which is going to become zero as we approach correct solution
"""
return y2-y2_old-dt*((m*y1)-(n*y2*y3)-(p*((y2)**2)))
def f3(y3_old,y1,y2,y3,dt):
"""
Function takes values of y1,y3_old,y2,y3 and dt and returning value which is going to become zero as we approach correct solution
"""
return y3-y3_old-dt*(p*((y2)**2))
def jacobian(y1_old,y2_old,y3_old,y1,y2,y3,dt):
"""
Defining Jacobian matrix
"""
J=np.zeros((3,3)) #2D Array defined
h=1e-6 #Very small value chosen to calculate Jacobian matrix values
"""
Since,Jacobian matrix values is calculated Numerically as
J=(f(x+h)-f(x))/h
Implimenting same below
"""
#Row 1
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
#Defining Initial guess value
y1_old=1
y2_old=0
y3_old=0
#Defining numpy column vector for y_old
y_old=np.zeros((3,1))
#Defining numpy column vector for y_old
y_guess=np.zeros((3,1))
#Defining numpy column vector for F
F=np.copy(y_old) #This means,F value will be same as y_old initially,but change in value of y_old will not affect values in F
total_t=600 #Total time 10 min=600 sec
dt=0.01 #Step size chosen
t=np.arange(0,total_t,dt) #Array of t is defined
Y=np.zeros((len(t),3)) #Array of final solution after 10 min
timestep_iter=0
#Starting Time marching loop
for i in range(0,len(t)):
y_old[0]=y1_old
y_old[1]=y2_old
y_old[2]=y3_old
"""
Initially value at previous iteration is chosen as y_guess and is equal to y_old
"""
y_guess[0]=y1_old
y_guess[1]=y2_old
y_guess[2]=y3_old
#Starting the newton raphson loop
iter=0
#Defining error
error=1e9
#Parameters for Newton raphson solver
alpha=1 #Relaxation parameter
tol=1e-12 #Tolerance value
"""
While for getting solution value till error>tol value
"""
while error>tol :
#For computing Jacobian y_old,y_guess values are used
J=jacobian(y_old[0],y_old[1],y_old[2],y_guess[0],y_guess[1],y_guess[2],dt)
#For computing F vector y_old,y_guess values are used
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)
#Computing new value y_new at next iteration
y_new=y_guess-(alpha*np.matmul(inv(J),F))
#Computing max error absolute value
error=np.max(np.abs(y_new-y_guess))
#Updating the new values
y_guess=y_new
#Log message
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
#Updating value of present time step to previous values
y1_old=y_new[0]
y2_old=y_new[1]
y3_old=y_new[2]
"""
Final solution after every time step getting stored in Y
"""
Y[i][0]=y_new[0]
Y[i][1]=y_new[1]
Y[i][2]=y_new[2]
log_message='Time step Iteration = {0} y1_end={1} y2_end={2} y3_end={3}'.format(timestep_iter,y_new[0],y_new[1],y_new[2])
print(log_message)
timestep_iter=timestep_iter+1
#print(Y)
plt.plot(t,Y[:,0],label='y1')
plt.plot(t,Y[:,1],label='y2')
plt.plot(t,Y[:,2],label='y3')
plt.xlabel('Time in sec ')
plt.ylabel('Change in values of y1,y2,y3')
plt.legend(loc="upper right")
plt.title('Plot showing trend of change in value of y1,y2,y3 w.r.t Time ')
plt.show()#Reason for choosing particular time step size:
As it is a Implicit method,so with change in step size solution will not become unstable(this method is stable),but only difference being smaller the step size less will be the difference between Numerical solution and Analytical Solution.
We will get more accurate solution only.As solution need to be found out for 10 minutes,so choosing very small step size will take longer time.Therefore choosing optimal step size that can give us benefit in accuracy and simulation time should be used.