Week 6 - Multivariate Newton Rhapson Solver
AIM :- TO SOLVE SET OF ORDINARY DIFFERENTIAL EQUATION USING MULTIVARIATE NETWON-RAPHSON TECHNIQUE.
INTRODUCTION :-
In this project, I have basically solved a set of ordinary differential equation using multivariate Newton-Raphson technique. Here, My aim was to analyse the effect of time- step on the solution of set of equation of ordinary differential equation for different variable. So, keeping in mind this concept, here I have made a program to find the solution of set of O.D.E equations containing three variable and plot the effect of time-step over the solution of this variables. This type of equations have wide range of application in engineering world where it is used in chemical kinetics in order to calculate the rate of reaction for formation of products where it includes enthalphy, entropy related with the reactions during combustion process. Apart from it, it has applications in heat conduction analysis like in design of heat spreaders in microelectronics, In combined heat conduction and convection like design of heating and cooling chambers.
THEORY :-
In order to solve the set of coupled odinary differential equation, We require a special technique which can solve such types of equations because this O.D.E are stiff equations which is very tough to solve.
Hence we either require a special method or special solvers provided by the software like MATLAB whcih can solve such equations. Here, in this project I have work on solver which can solve such type of coupled equations. Since python is very user-friendly programming language and have a fast respone output characteristics, I have chosen it to built a solver. For this ,I have used mulltivariate Newton-Raphson method tehnique which can solve the equations.
Mutivariate Newton-Raphson method :-
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`(x1, x2, x3, x4..........xn-1, xn)=0
`f_2`(x1, x2, x3, x4..........xn-1, xn)=0
`f_3`(x1, x2, x3, x4..........xn-1, xn)=0
...
`f_(n-1)`(x1, x2, x3, x4..........xn-1, xn)=0
`f_n`(x1, x2, x3, x4..........xn-1, xn)=0
The formula for Multivariate Newton-Raphson method is:-
`Y_(n ew)= Y_(guess) - alpha*J^(-1)F`
where 'J' stands for Jacobian matrix and it is in the following way :-
J= `[[(delf_1)/(delx_1), (delf_1)/(delx_2) (delf_1)/(delx_3) ...(delf_1)/(delx_n)], [(delf_2)/(delx_1), (delf_2)/(delx_2) (delf_2)/(delx_3)...(delf_2)/(delx_n)], [(delf_n)/(delx_1),(delf_n)/(delx_2) (delf_n)/(delx_3)... (delf_n)/(delx_n)]]`
and, `alpha` is known as relaxation factor. `Y_(n ew)` is a row matrix of n variables.
Here in the Program, we have considered following given ordinary differential equations :-
`(dy_1)/(dt)= -0.04y_1 + 10^4y_2y_3.`
`(dy_2)/dt=0.04y_1 - 10^4y_2y_3 - 3**10^7y_2^2.`
`dy_3/(dt)=3**10^7y_2^2.`
The above equations have been formulated using Euler implict/backward differencing method in the following way:-
`y_1^(n)-y_1^(n-1) - Deltat(-0.04**y_1^n + 10^4**y_2^n**y_3^n)=0.`
`y_2^(n)-y_1^(n-1) - Deltat(0.04**y_1^n - 10^4**y_2^n**y_3^n -3**10^7**y_2^n)=0.`
`y_3^n-y_2^(n-1) -Deltat(3**10^7y_2^n)=0.`
The given O.D.E have been formulated in above form so that it could be make convenient for programming methodology.
"""
PROGRAM TO SOLVE ORDINARY DIFFERENTIAL EQUATION USING MULTIVARIATE NEWTON-RAPHSON TECHINQUE.
HERE IN THIS PROGRAM, WE ARE SOLVING SET OF COUPLED ORDINARY DIFFERENTIAL EQUATION WHICH IS GIVEN BELOW :-
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
"""
import matplotlib.pyplot as plt
import numpy as np
from numpy.linalg import inv
def f1(y1_old, y1, y2, y3, dt): # function defined for first ordinary differential equation.
return y1_old + dt*(-(0.04)*y1 + (10^4)*y2*y3) - y1
def f2(y2_old, y1, y2, y3, dt): # function defined for second ordinary differential equation.
return y2_old + dt*((0.04)*y1 - (10^4)*y2*y3 - 3*(10^7)*(pow(y2,2))) - y2
def f3(y3_old, y1, y2, y3, dt): # function defined for third ordinary differential equation.
return y3_old + dt*(3*(10^7)*(pow(y2,2))) - y3
def Jacobian(y_old, y1, y2, y3, dt): # function defined for Jacobian matrix of order 3x3.
"""
Defining the jacobian.
"""
d=1e-4 # defined grid size of value 0.1
J= np.ones((3,3)) # intialize Jacobian matrix of order 3x3 having value 1.
# Row1
J[0,0]= (f1(y_old[0], y1+d, y2, y3, dt) - f1(y_old[0], y1, y2, y3, dt))/d;
J[0,1]= (f1(y_old[1], y1, y2+d, y3, dt) - f1(y_old[1], y1, y2, y3, dt))/d;
J[0,2]= (f1(y_old[2], y1, y2, y3+d, dt) - f1(y_old[2], y1, y2, y3, dt))/d;
#row2
J[1,0]= (f2(y_old[0], y1+d, y2, y3, dt) - f2(y_old[0], y1, y2, y3, dt))/d;
J[1,1]= (f2(y_old[1], y1, y2+d, y3, dt) - f2(y_old[1], y1, y2, y3, dt))/d;
J[1,2]= (f2(y_old[2], y1, y2, y3+d, dt) - f2(y_old[2], y1, y2, y3, dt))/d;
#row3
J[2,0]= (f3(y_old[0], y1+d, y2, y3, dt) - f3(y_old[0], y1, y2, y3, dt))/d;
J[2,1]= (f3(y_old[1], y1, y2+d, y3, dt) - f3(y_old[1], y1, y2, y3, dt))/d;
J[2,2]= (f3(y_old[2], y1, y2, y3+d, dt) - f3(y_old[2], y1, y2, y3, dt))/d;
return J
# defining the intial guess values of y1,y2,y3
y1= 1
y2= 0
y3= 0
# intial error guess value
err=9e9
# define a matrix of order 3x1 and assigning them value of variable y1, y2, y3.
y_old=np.ones((3,1))
y_old[0]=y1;
y_old[1]=y2;
y_old[2]=y3;
# define a matrix of order 3x1 and assigning them value of 0.
y_new=np.zeros((3,1))
y_new[0]=1
alpha=1.2 # relaxation-factor.
tol=1e-15 # tolerance value.
time_steps=60001 # total time-steps.
dt=0.01 # time-step value.
t=np.arange(0,time_steps,1)
# here defined matrix of order of length of 't' for three variables given below to store value of y after each time integration.
y1_upd=np.zeros(time_steps)
y2_upd=np.zeros(time_steps)
y3_upd=np.zeros(time_steps)
iter=1 # iteration counter to count number of iterations.
# starting newton-rhapson method
for i in range(1,time_steps): # time - integration loop.
y_prev=np.ones((3,1))
y_prev[1]=0
y_prev[2]=0
F=np.copy(y_prev)
while (err>tol): # convergence loop.
J=Jacobian(y_old, y_prev[0], y_prev[1], y_prev[2], dt)
F[0]=f1(y_old[0], y_prev[0], y_prev[1], y_prev[2], dt)
F[1]=f2(y_old[1], y_prev[0], y_prev[1], y_prev[2], dt)
F[2]=f3(y_old[2], y_prev[0], y_prev[1], y_prev[2], dt)
y_new= y_prev - alpha*(np.matmul(inv(J),F))
err= np.max(np.abs(y_new - y_prev))
y_prev=y_new
y_old=y_new
y1_upd[i]=(y_old[0])
y2_upd[i]=(y_old[1])
y3_upd[i]=(y_old[2])
iter=iter+1
err=1 # assignment of variable after each loop of convergence.
log_message='iteration no = {0} y1={1} y2={2} y2={3}' .format(iter, y_new[0], y_new[1], y_new[2])
print(log_message)
plt.plot(t,y1_upd,'b',label='VALUE OF Y1')
plt.plot(t,y2_upd,'r',label='VALUE OF Y2')
plt.plot(t,y3_upd,'g',label='VALUE OF Y3')
plt.legend(loc="upper right")
plt.xlabel('NUMBER OF TIME-STEPS(TIME-STEP VALUE=0.01)')
plt.ylabel('VALUE OF DIFFERENT VARIABLE OF ORDINARY DIFFERENTIAL EQUATION')
plt.title('SOLUTION USING MULTIVARIATE NEWTON RAPHSON METHOD WITH TIME INTEGRATION',weight='bold')
plt.show()
Here, while observing the different output plot, I observed that the best output which I obtained was at value of time-step 0.01 . Hence I concluded that the best time-step for getting the right output result is 0.01.
This output plot was generated after iterations of 60000 time steps. So, the output result for value of y1, y2 and y3 after 60000 iterations are -
Note :- The simulation have been computed for approx 10 minutes.