Week 6 - Multivariate Newton Rhapson Solver
Objective: To solve a given set of Ordinary Differential equations using the Multi-Variate Newton Raphson Method.
Given: The set of ODE's are given below:
`dy_1/dt= - 0.04*y_1+10^4*y2*y3`
`dy_2/dt=0.04*y_1-10^4*y_2*y_3-3*10^7*y_2^2`
`dy_3/dt=3*10^7*y_2^2`
The jacobian should be estimated numerically and not analytically. The initial conditions are 1 for y1 and 0 for y2 and y3 respectively.
Take alpha as 1. Run the simulation for 10 minutes.
Theory:
`"Multi-Variate Newton Raphson Method"`
As we know, the Multivariate Newton-Raphson method is a direct extension of the single variable Newton-Raphson method. Where the single variable Newton-Raphson method solved
`f(x)=0`, the multivariate version will solve a system of n equations of the form:
`f_1(x_1,x_2,x_3.....,x_(n-1),x_n) = 0`
`f_2(x_1,x_2,x_3.....,x_(n-1),x_n) = 0`
`f_3(x_1,x_2,x_3.....,x_(n-1),x_n) = 0`
.
.
.
.
`f_(n-1)(x_1,x_2,x_3.....,x_(n-1),x_n) = 0`
`f_n(x_1,x_2,x_3.....,x_(n-1),x_n) = 0`
We will adopt the shorthand notation for the equation.
`f(x)=0`
According to Multi-variate Newton Raphson Method, to get the solution out of coupled non-linear ODE below equation we use:
`y_("new")=y_("old")-alpha*J^(-1)*F`
where, alpha - Relaxation factor
Vector F will be written as,
`f=[[f_1],[f_2],[f_3],[.],[.],[f_n]]`
Old value or Guess value will be written as,
`y_("old")=[[y_("1 old")],[y_("2 old")],[y_("3 old")],[.],[.],[y_("n old")]]`
Jacobian Matrix,
`J=[[(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_2)/dy_3],[.],[.]]`
where, nos. of column= Nos. of independent variables
nos. of Rows= Nos. of Functions
Solution:
Python Code:
"""
Program to solve non-linear coupled ODE with multi-variate Newton Raphson Solver
System of coupled non-linear ODE is:
-0.04 * Y1 + 10^4 * Y2*Y3 = 0
0.04 * Y1 - 10^4 *Y2*Y3 - 3 * 10^7* Y2^2 = 0
3 * 10^7* Y2^2 = 0
By- Gaurav V Kharwade || Skill-Lync:2020
"""
import numpy.matlib
import numpy as np
import numpy.linalg
import matplotlib.pyplot as plt
def f1(y1,y2,y3,y1_old,dt): # First non-linear equation
return (y1_old + dt*(-0.04*y1 + 1e4*y2*y3) - y1)
def f2(y1,y2,y3,y2_old,dt): # Second non-linear equation
return (y2_old + dt*(0.04*y1 - (1e4*y2*y3) - (3*1e7*y2*y2)) - y2)
def f3(y1,y2,y3,y3_old,dt): # Third non-linear equation
return (y3_old + dt*(3e7*y2*y2) - y3)
def jacobian(y1,y2,y3,y1_old,y2_old,y3_old,dt): # To define Jacobian matrix
#To define 3 * 3 matrix creation
j= np.ones((3,3))
# Step-size
h = 1e-6
"""
For estimation of jacobian matrix we have used numerical method(BDS) instead analytical method
"""
# Row 1
j[0,0]= (f1(y1+h,y2,y3,y1_old,dt)-f1(y1,y2,y3,y1_old,dt))/h
j[0,1]= (f1(y1,y2+h,y3,y1_old,dt)-f1(y1,y2,y3,y1_old,dt))/h
j[0,2]= (f1(y1,y2,y3+h,y1_old,dt)-f1(y1,y2,y3,y1_old,dt))/h
# Row 2
j[1,0]= (f2(y1+h,y2,y3,y2_old,dt)-f2(y1,y2,y3,y2_old,dt))/h
j[1,1]= (f2(y1,y2+h,y3,y2_old,dt)-f2(y1,y2,y3,y2_old,dt))/h
j[1,2]= (f2(y1,y2,y3+h,y2_old,dt)-f2(y1,y2,y3,y2_old,dt))/h
# Row 3
j[2,0]= (f3(y1+h,y2,y3,y3_old,dt)-f3(y1,y2,y3,y3_old,dt))/h
j[2,1]= (f3(y1,y2+h,y3,y3_old,dt)-f3(y1,y2,y3,y3_old,dt))/h
j[2,2]= (f3(y1,y2,y3+h,y3_old,dt)-f3(y1,y2,y3,y3_old,dt))/h
return j
# Initialisation
y1_old= 1
y2_old= 0
y3_old= 0
# Old value assignment
y_old= np.ones((3,1))
y_old[0]= y1_old
y_old[1]= y2_old
y_old[2]= y3_old
# Column vector of function F
F= np.copy(y_old)
# Relaxation factor
alpha= 1
# Total time, T=10 mins, t= 10 * 60
t= 600
# Time step
dt= 0.1
n= int(t/dt)
# Array of updated values
y_new1= []
y_new2= []
y_new3= []
it= []
error_updated= []
# Initialising iteration
iter= 1
# For loop, to time stepping
for i in range(0,n):
# Newton Raphson criteria
tol= 9e-6
error= 9e9
# Newton Raphson Solver
while error > tol:
# Jacobi Matrix
j= jacobian(y_old[0],y_old[1],y_old[2],y1_old,y2_old,y3_old,dt)
# F column vector
F[0]= f1(y_old[0],y_old[1],y_old[2],y1_old,dt);
F[1]= f2(y_old[0],y_old[1],y_old[2],y2_old,dt);
F[2]= f3(y_old[0],y_old[1],y_old[2],y3_old,dt);
# Equation for solver
y_new = y_old - alpha * np.matmul(numpy.linalg.inv(j), F);
# Error finding
error = np.max(np.abs(y_new - y_old));
print(error)
# Old value updation
y_old = y_new
# To print on console
log= 'iter= {0} error= {1} y1= {2} y2= {3} y3= {4} '.format(iter,error, y_new[0], y_new[1], y_new[2])
print(log)
y_old = y_new;
y_new1.append(y_new[0])
y_new2.append(y_new[1])
y_new3.append(y_new[2])
error_updated.append(error)
iter= iter + 1
it.append(iter)
# Old value updation
y1_old= y_new[0]
y2_old= y_new[1]
y3_old= y_new[2]
# To create a plot of error dropping
plt.plot(it,error_updated,'--',color='red')
plt.yscale("log")
plt.title('Error dropping characteristics')
plt.xlabel('Iteration nos.')
plt.ylabel('Error value')
plt.grid('on')
# To create a figure with 3 rowa and 1 column of subplot
fig, ax= plt.subplots(3)
# For convergence of y1
ax[0].plot(it,y_new1)
ax[0].set_title('Convergence of y1')
ax[0].grid('on')
# For convergence of y1
ax[1].semilogy(it,y_new2)
ax[1].set_title('Convergence of y2')
ax[1].grid('on')
# For convergence of y3
ax[2].plot(it,y_new3)
ax[2].set_title('Convergence of y3')
ax[2].grid('on')
fig.tight_layout()
plt.show()
Results:
After running the program the results we get is as below:
We can see above the error values that we have targeted have achieved and the solution for this ODE has shown above.
y1= 0.39997579, y2= 2.63177e-6 & y3= 0.600021
We can see below, at each and every iterations we are getting closer to an accurate value.
In this plot, we can see how the errors are getting dropped at every iteration in order to achieve a solution with high accuracy.
I ran the program by using a time step much lower than 0.1, but it has been observed that there was no major variation in the solution obtained. Hence, choosing 0.1 as time step will be a better choice as it helps us to get solutions in less time and with high accuracy also.