Week 6 - Multivariate Newton Rhapson Solver
"""
program to solve the system of
coupled non-linear equation
dy_1/dt = -0.04y_1 + 10^4y_2y_3
dy_1/dt = 0.04y_1 - 10^4y_2y_3 - 3.10^7y_2^2
dy_3/dt = 3.10^7y_2^2
Modulus required
1. numpy is required for matrix multiplication
2. inv is required for computing the inverse of a matrix
"""
import numpy as np
from numpy.linalg import inv
import matplotlib.pyplot as plt
# For 1st equation
def f1(y1_old, y1, y2, y3, dt):
#return y1_old + dt*(-0.04*y1 + (pow(10,4)))*y2*y3 - y1
return y1_old + dt*(-0.04*y1 + 10**4*y2*y3) - y1
# For 2nd equation
def f2(y2_old, y1, y2, y3, dt):
#return y2_old + dt*(0.04*y1 - (pow*(10,4))*y2*y3 - 3*(pow(10,7))*(pow(y2,2)) - y2
return y2_old + (0.04*y1 - 10**4*y2*y3 - 3*10**7*y2**2)*dt - y2
# For 3rd equation
def f3(y3_old, y1, y2, y3, dt):
#return y3_old + dt*(3*(pow(10,7))*(pow(y2,2)) - y3
return y3_old + dt*(3*10**7*y2**2) - y3
def jacobian(y1_old, y2_old, y3_old, y1, y2, y3, dt):
J = np.ones((3,3))
dy = 1e-4
J[0,0] = (f1(y1_old, y1+dy, y2, y3, dt) - f1(y1_old, y1, y2, y3, dt))/dy
J[0,1] = (f1(y1_old, y1, y2+dy, y3, dt) - f1(y1_old, y1, y2, y3, dt))/dy
J[0,2] = (f1(y1_old, y1, y2, y3+dy, dt) - f1(y1_old, y1, y2, y3, dt))/dy
J[1,0] = (f2(y2_old, y1+dy, y2, y3, dt) - f2(y2_old, y1, y2, y3, dt))/dy
J[1,1] = (f2(y2_old, y1, y2+dy, y3, dt) - f2(y2_old, y1, y2, y3, dt))/dy
J[1,2] = (f2(y2_old, y1, y2, y3+dy, dt) - f2(y2_old, y1, y2, y3, dt))/dy
J[2,0] = (f3(y3_old, y1+dy, y2, y3, dt) - f3(y3_old, y1, y2, y3, dt))/dy
J[2,1] = (f3(y3_old, y1, y2+dy, y3, dt) - f3(y3_old, y1, y2, y3, dt))/dy
J[2,2] = (f3(y3_old, y1, y2, y3+dy, dt) - f3(y3_old, y1, y2, y3, dt))/dy
return J
# Defining the initial guess values
y1 = 1
y2 = 0
y3 = 0
# numpy column vector
y_old = np.ones((3,1))
y_old[0] = y1
y_old[1] = y2
y_old[2] = y3
# numpy column vector for F
# why not F = X_old?
F = np.copy(y_old)
# parameters for newton rhapson solver
alpha = 1
dt = 0.1
t = 600
n = int(t/dt)
iter = 1
y_new = np.ones((3,1)) # New Value
# For plot, I have store the values of y1,y2, y3 and iteration in array format using append
y_1 = []
y_2 = []
y_3 = []
iter_array = []
# starting the Newton Rhapson loop
for i in range(0,n):
error = 9e9
tol = 1e-6
while error>tol:
# Computing the jacobian - old values are used
J = jacobian(y_old[0], y_old[1], y_old[2], y1, y2, y3, dt)
# Computing the F vector - old values are used
F[0] = f1(y1, y_old[0], y_old[1], y_old[2], dt)
F[1] = f2(y2, y_old[0], y_old[1], y_old[2], dt)
F[2] = f3(y3, y_old[0], y_old[1], y_old[2], dt)
y_new = y_old - alpha*(np.matmul(inv(J),F))
error = np.max(np.abs(y_new - y_old))
#print(error)
y_old = y_new
# log message
log_message = 'iteration={0} y1={1} y2={2} y3={3}'.format(iter, y_old[0], y_old[1], y_old[2])
print(log_message)
iter = iter+1
# update values (For plots - stored the values)
y_1.append(y_new[0])
y_2.append(y_new[1])
y_3.append(y_new[2])
iter_array.append(iter)
# final values
y1 = y_new[0]
y2 = y_new[1]
y3 = y_new[2]
# Plots
plt.plot(iter_array,y_1, '--',linewidth = 1.5,color = 'purple')
plt.plot(iter_array,y_2,'-',linewidth = 1.5,color = 'red')
plt.plot(iter_array,y_3,linewidth = 1.5,color = 'yellow')
plt.legend()
plt.title('Y values for dt = 0.1')
plt.xlabel('Iteration number')
plt.ylabel('Y values')
plt.grid()
plt.show()AIM - To solve the set of equations with multiple variables using the newton Rhapson solver.
INTRODUCTION
- While solving a combustion mechanism,the rates of production/depletion of different species are calculated in the form of ordinary differential equations (ODE's)
- The rate of production/dissipation of a particular species depends on the concentration of the species itself as well as concentration of other species. such set of ODE's is called a system of coupled differential equations.
- They cannot be solved directly through integration;hence numerical methods are used to integrate and solve the system of coupled differential equations. One such method is the multivariate Newton-Raphson method, which is an extention of the univariate Newton-Raphson method.
- The Newton-Raphson method solves the root of a function through iterative techniques. For a univariate system, we can find the root of the function through guess values as:
`y_n = y_{guess} - alpha [ f_(y_{guess}) / f_(y_{guess})]`
Key differences between Implict and Explicit methods
Implicit method will be stable for larger values of step size (but not accurate)
Explicit methods will become un-stable if time step increases above a critical method.
For Stiff systems implicit method using Backward Difference Method is preferred.
Given data -
`dy_1/dt = -0.04y_1 + 10^4y_2y_3`\
`dy_1/dt = 0.04y_1 - 10^4y_2y_3 - 3.10^7y_2^2`\
`dy_3/dt = 3.10^7y_2^2`
Here i am going to putting all the species together i actually have a ordinary differential equation for simplicity and also i am going to know down 3 ordinary differential equation dy1/dt , dy2/dt, dy3/dt
i am solving for y1,y2 and y3 or we can think in this way that there are species it may be y1 can be carbon dioxide, y2 can be mono-oxide and so on. I am basically saying that how the concentration of y1,y1 and y3 are going to change with time.
Now, if we look at attention to this ODE system we can called as coupled ODE system. if we take a look at 1st ODE we are solving for y1 and we can see in the right hand side we have y2 and y3 are involved but we have seprate ODE for y2 , which means we can not solve the 1st ODE seprately and use those results to solve 2nd ODE.
The general layout of the code is given below :
- Outer time integration loop
- Assembling the ODE system
- Newton Rhapson solver till convergence
- Get new values and update them
- End of the outer loop
Newton Rhapson method is implementedat every time step for the given simulation time.
simulatiion time = 10 mintues
alpha = 1
The initial conditions given are,
y1 = 1
y2 = 0
y3 = 0
MATHEMATICAL FORMULATION : -
1.INTEGRATION USING IMPLICIT /BACKWARD DIFFERENCING METHOD :
Let us integrate the first equation of the coupled non linear system making use of the implicit /backward differencing method and also rearrange the equation:
Now, all the equations are discretized using the backward differencing as,
`dy_1/dt = -0.04y_1^n + 10^4 (y_2 y_3)^n`
`(y_1^n - y_1^{n-1})/{dt} = -0.04y_1^n + 10^4(y_2y_3)^n`
`y_1^n =y_1^{n-1} + dt ( -0.04y_1^n + 10^4(y_2y_3)^n)`
`y_1^{n-1} + dt ( -0.04y_1^n + 10^4(y_2y_3)^n) - y_1^n = 0`
2nd equation be like
`dy_2/dt = 0.04y_1 - 10^4y_2y_3 - 3*10^7y_2^2`
`(y_2^n - y_2^{n-1})/{dt} = 0.04y_1 - 10^4y_2y_3 - 3*10^7y_2^2``
`y_2^n = y_2^{n-1} + dt*( 0.04y_1 - 10^4y_2y_3 - 3*10^7y_2^2)`
`y_2^{n-1} + dt*( 0.04y_1 - 10^4y_2y_3 - 3*10^7y_2^2)- y_2^n = 0`
Meanwhile 3rd equation be like,
`dy_3/dt = 3.10^7y_2^2`
`(y_3^n - y_3^{n-1})/{dt} = 3.10^7y_2^2`
`y_3^n =y_3^{n-1}dt*( 3.10^7y_2^2)`
`y_3^{n-1}dt*( 3.10^7y_2^2)-y_3^n = 0`
2. SETTING UP A COLUMN MATRIX
The first equation when we solve relates to a vector and we get a set of vectors for each and every time step. The function is evaluated for the current time step. For a time step n = 1, we get the following vector
`f_1 = y_1^0 + dt ( -0.04y_1^1 + 10^4(y_2y_3)^1) - y_1^1 = 0`
`f_2 = y_2^0+ dt*( 0.04y_1 - 10^4y_2y_3 - 3*10^7y_2^2)- y_2^1 = 0``
`f_3 = y_3^3 + dt*( 3.10^7y_2^2)-y_3^1 = 0``
All those elements are now packaged into a column vector F given by :
F = `F = [[f_1],[f_2],[f_3]]`
3. MULTIVARIATE NEWTON RAPHSON TECHNIQUE :
We have multiple equations and we we make use of multivariate newton raphson technique. We start with an initial array of values and we put them in a column matrix. The formula is given below :
`y^{New} = y^{old} - alpha J^-1 F`
`alpha` = Relaxation Factor
F = function (column vector)
J = jacobian Matrix
Representing all the values as a column matrix we get the following :
`[[y_1^{1,New}],[y_2^{2,New}],[y_3^{3,New}]]` = `[[y_1^{1,old}],[y_2^{1,old}],[y_3^{1,old}]]` - `alpha``[[{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_1],[f_2],[f_3]]`
All possible derivatives for the matrix are written making use of the Jacobian
For y initially we start with guess values and then get new value for y.
The new values of y are again stored as initial guess values and the process continues to the next iteration and so on.
4. NUMERICAL DIFFERENTIATION OF THE JACOBIAN MATRIX :
For differentiating our mathematics function numerically inside the Jacobian Matrix we use forward difference formula. For a function f1 the formula is provided as show below:
`{df_1}/{dy_1} ={f(y1+dy) - f(y1)}/{dy}`
Similarly we differentiate for all possible derivatives within the Jacobian Matrix. The forward differencing is only possible for a small incremental value of dy. The value must be made as small as possible to provide better results
Coding :
For solving the actual problem using Python, we proceed as usual and solve the quadratic equations first by defining functions for each equations. The jacobian matrix is then differentiated numerically,by discretizing using for forward differencing scheme.
Dummy values are assumed ina array which will be used in the time loop required for the problem . In our case, we are asked to simulate th system for 10 mins. i.e, for 600 seconds. So the solution arrays defined are given lengths according to the total time step required to solve the system.
The error and initial values are defined before the time loop starts and are accordingly updated at the end of the loop for the next iteration to take place.
Conclusion -
For the time step dt, we have assume the best value that gives a fairly accurate results for the problem. if we take dt as 0.5 there are 1200 steps in the system and hence results are sharp and not that accurate. But if we give a smaller dt as 0.1, we have too many iterations to complete in an already complex multivariate solver, eventhough the solution result comes out smooth. Hence an average dt value is assumed in this case to give an assumed in the case to give an accurate representation of data an also sufficient time efficient.
The results show the variations in values and their domains getting smaller due to increases in time step size . Hence the ideal time step for this case is aasumed to be 0.1
More iterations needs to be run for better accuracy in the results for simulation. Having the value of dt<0.1 provides more accurate results for our simulation but there is a necessity for bringing down the time required for computation.
The unconditional stability of the solution is due to the fact that we are using an implicit scheme here which is stable irrespective of the timestep value. where as the stability of the explicit scheme depends on the value of time step.
Through explicit schemes are easy to solve compared to the implicit schemes.
From the Given attached graphs we can see that the approximate values of y1,y2,y3 are changing continuously with respect to time. For the values of y3 we cannot observe the change because the value of y2 are really small .The values changed rapidly as the solution marches forward in time. This behaviour shows the stiff characteristics of the system.
The Value of timestep we used was 0.1 but changing the value of the timestep will not affect the results very much. The reason for this we are using a implict scheme. This only issue we will be facing while changing the time step would be the computational time.
OUTPUT :-
1. As we move forward in time the value of function varies rapidly.
2. Change in the time step doesn't make any big difference in final solution.
3. Small value of time step taken longer time (iterations) to converged results which may increase cost.
4. considering small timestep increases number of nodes which gives more precise solution curve.
5. As implicit scheme is considered we get stable solution but setting implicit scheme is more hectic and more complex.