Multivariate Newton Rhapson Technique
Objective-To solve the set of ODE equations using Multivariate Newton Rhapson Technique.
The given set of ODE's are:
`(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`
The problem is discretized using Implicit Euler Method/Backward Differencing.
`(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 above equations in tems of RHS-LHS=0 we get the following non-linear equations,
`f1=y1^(n-1)+Deltat*(-0.04*y1^n+10^4*y2^n*y3^n)-y1^n`
`f2=y2^(n-1)+Deltat(0.04*y1-10^4*y2*y3-3*10^7*y2^2)^n-y2^n`
`f3=y3^(n-1)+Deltat(3*10^7*y2^2)^n-y3^n`
`y^n`=The value at current time step
`y^(n-1)`The value at previous time step
The Newton Rhapson formula is given as
`Y_(New)=Y_(guess)-alpha*J^-1.F`
The initial values is taken as:
`Y1_(guess)=1 ,Y2_(guess)=0,Y3_(guess)=0`
Matrix F is the column matrix of the functions,
`F=[(f1),(f2),(f3)]`
Jacobian matrix J is given as,
`J=[((delf1)/(dely1),(delf1)/(dely2),(delf1)/(dely3)),((delf2)/(dely1),(delf2)/(dely2),(delf2)/(dely3)),((delf3)/(dely1),(delf3)/(dely2),(delf3)/(dely3))]`
`alpha=1`
Coding approach:
1. The derivatives in jacobian matrix is solved using numerical technique,
`(delf)/(dely)=(f(y+h)-f(y))/h`
2. The initial guess values are taken and the ODE system solved till convergence using a while loop.
3. The new values are updated for each time level using a for loop.The simulation is run for 10 min i.e 600 seconds using a time step of 0.5
Program Structure:
Outer time integration loop
Assembling the ODE system
Run the Newton Rhapson till convergence
Get new values+update
End time loop
Python Code:
import numpy as np
import matplotlib.pyplot as plt
from numpy.linalg import inv
def f1(y1_old,y1,y2,y3,dt):
return y1_old + dt*(-0.04*y1+1e4*y2*y3)-y1
def f2(y2_old,y1,y2,y3,dt):
return y2_old + dt*(0.04*y1-1e4*y2*y3-3e7*y2**2)-y2
def f3(y3_old,y1,y2,y3,dt):
return y3_old + dt*(3e7*y2**2)-y3
def jacobian(y1_old,y2_old,y3_old,y1,y2,y3,dt):
J=np.ones((3,3))
h=1e-7
#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
#row 2
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
#row 3
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 values
y1=1
y2=0
y3=0
yguess=np.ones((3,1))
yguess[0]=y1
yguess[1]=y2
yguess[2]=y3
y_old=np.ones((3,1))
y_old[0]=y1
y_old[1]=y2
y_old[2]=y3
F=np.copy(yguess)
alpha=1
tol=1e-6
dt=0.5
y_final_1=[]
y_final_2=[]
y_final_3=[]
time=np.arange(0,600,dt)
iter=1
for t in np.arange(0,600,dt):
error=9e9
while error>tol:
J=jacobian(y_old[0],y_old[1],y_old[2],yguess[0],yguess[1],yguess[2],dt)
F[0]=f1(y_old[0],yguess[0],yguess[1],yguess[2],dt)
F[1]=f2(y_old[1],yguess[0],yguess[1],yguess[2],dt)
F[2]=f3(y_old[2],yguess[0],yguess[1],yguess[2],dt)
ynew=yguess - alpha*np.matmul(inv(J),F)
error=np.max(np.abs(ynew-yguess))
yguess=ynew
iter=iter+1
log_message='iteration #={0} y1={1} y2={2} y3={3}' .format(iter,ynew[0],ynew[1],ynew[2])
print(log_message)
y_old=ynew
y_final_1.append(ynew[0])
y_final_2.append(ynew[1])
y_final_3.append(ynew[2])
plt.plot(time,y_final_1,label='y1')
plt.plot(time,y_final_2,label='y2')
plt.plot(time,y_final_3,label='y3')
plt.xlabel('time in sec')
plt.ylabel('Y')
plt.legend()
plt.grid()
plt.show() Results:
The below plot show the variation y with respect to time.A time step choosen for this is 0.5s and simulation was run for 10 minutes.
Time step 0.5 s
Time step 1 s
Time step 30 s
From the above results it can be seen that as we decrease the time step, the number of iterations increases therby increasing the compution time.There is not much difference in results, but with larger time step the accuracy has been compromised.
A optimum time step of 0.5 s was choosen as decreasing the time step further from 0.5 s did not change the solution much.
Conclusion:
1. It is seen that when the time step is varied the solution is stable throughout and harldy there is any change in the solution.The reason for this is the use of implicit method.Imcplicit methods are unconditionally stable.
2. So with large time time step we could get the solution,which is not possible by using explicit method.But the accuracy is decreased with larger time step.