Newton Rhapson Solver for Coupled Non-linear partial differential equations

Objective:
In the project, the aim is to script a Multivariate Newton Rhapson Solver in Python for solving the coupled non-linear partial differential equation. The set equations used for solving this particular case is shown below:

`(dy_1)/dt = -0.04y_1 + 10^4 y_2 y_3`

`(dy_2)/dt = 0.04y_1 - 10^4 y_2 y_3 - 3*10^7 y_2^2`

`(dy_3)/dt = 3*10^7y_2^2`

Solution:
When observing the given system of equations, we can say that they are coupled and non-linear. We establish that they are coupled because, y1, y2, y3 can only be solved by solving the three equations together. We say that it is non-linear because the unknown elements in the equations are getting multiplied by themselves or each other. Here, the Newton Rhapson method is implemented for each time step set for the given simulation time. The difference is that, in Multivariate, the approximations are calculated in a matrix form as there are more than one unknowns. The general principle for the Multivariate Newton Rhapson is given below:

`Xn+1 = Xn − α ⋅ J −1 ⋅ F`

where `X_n` is the old value,

`X_(n+1)` is the new value which we are guessing with respect to the old value,

`alpha` is the relaxation factor,

`J ^-1` is the inverse Jacobian matrix and

F is a matrix that contains the given system of equations.

The matrix representation of the above principle is given below:`

In the above figure, `alpha` is multiplied to the Jacobian inverse and the function matrix term.

In the current challenge, we have calculated the Jacobian matrix numerically as analytical derivatives are not always available. The principle of the numerical Jacobian is expressed below:

`(df_i)/dx_j = (f(x_j +h) - f(x_j)) / h`

In the Newton Rhapson solver, the relaxation factor is taken as 1.0. We are also using Implicit Euler or otherwise known as Backward Differencing method for the outer time integration loop for
this challenge.

Backward Differencing method is a finite approximation method which uses the function values at x and x − h where h is the step size. The general principle for this method is mentioned below:

`f'(x) = (f(x) - f(x-h)) / h`

Given the step size (h), note that this formula uses the values at and the point at the previous step. As it moves in the backward direction, it is called the backward difference method.

The given initial conditions for the unknowns are:

`y_1 = 1`

`y_2 = 0`

`y_3 = 0`

Here, we are initializing the ODE's in the form;

`f'(x) = x_(old value) + dt.f(x) - x _(solved value)`

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 Outer Loop

The code for the solver is shown below:

'''
Program to solve system of coupled non-linear equations using Multivariate Newton Rhapso
n Method
-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
'''
import numpy as np
from numpy.linalg import inv
import matplotlib.pyplot as plt

#Step size initialization
dt = 0.6
t_start = 0
t_end = 600 

# Initial Guess Values
y1 = 1
y2 = 0
y3 = 0

def f1(y1_old, y2_old, y3_old, y1, y2, y3, dt):
return y1_old + dt*(-0.04*y1 + pow(10,4)*y2*y3) - y1

def f2(y1_old, y2_old, y3_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

def f3(y1_old, y2_old, y3_old, y1, y2, y3, dt):
return y3_old + dt*(3*pow(10,7)*pow(y2,2)) - y3

def jac(y1_old,y2_old,y3_old,y1, y2, y3, dt):
# Defining the Jacobian Matrix
J = np.ones((3,3))
d = 1e-8
# Row 1
J[0,0] = (f1(y1_old, y2_old, y3_old, y1+d, y2, y3, dt) - f1(y1_old, y2_old, y3_old, y1, y2, y3, dt))/d
J[0,1] = (f1(y1_old, y2_old, y3_old, y1, y2+d, y3, dt) - f1(y1_old, y2_old, y3_old, y1, y2, y3, dt))/d
J[0,2] = (f1(y1_old, y2_old, y3_old, y1, y2, y3+d, dt) - f1(y1_old, y2_old, y3_old, y1, y2, y3, dt))/d

# Row 2
J[1,0] = (f2(y1_old, y2_old, y3_old, y1+d, y2, y3, dt) - f2(y1_old, y2_old, y3_old, y1, y2, y3, dt))/d
J[1,1] = (f2(y1_old, y2_old, y3_old, y1, y2+d, y3, dt) - f2(y1_old, y2_old, y3_old, y1, y2, y3, dt))/d
J[1,2] = (f2(y1_old, y2_old, y3_old, y1, y2, y3+d, dt) - f2(y1_old, y2_old, y3_old, y1, y2, y3, dt))/d

# Row 3
J[2,0] = (f3(y1_old, y2_old, y3_old, y1+d, y2, y3, dt) - f3(y1_old, y2_old, y3_old, y1, y2, y3, dt))/d
J[2,1] = (f3(y1_old, y2_old, y3_old, y1, y2+d, y3, dt) - f3(y1_old, y2_old, y3_old, y1, y2, y3, dt))/d
J[2,2] = (f3(y1_old, y2_old, y3_old, y1, y2, y3+d, dt) - f3(y1_old, y2_old, y3_old, y1, y2, y3, dt))/d

return J

def Newton_Rhapson(y1, y2, y3, y1_guess, y2_guess, y3_guess, dt):
# Old Guess Values Column vector
Y_old = np.ones((3,1))
Y_old[0] = y1_guess
Y_old[1] = y2_guess
Y_old[2] = y3_guess
F = np.ones((3,1))

#Parameters for Newton Raphson solver
error = 9e9 #Initializing error
alpha = 1
tol = 1e-9
iter = 1

while (error > tol):
J = jac(y1, y2, y3, Y_old[0],Y_old[1],Y_old[2], dt)

#Computing the main functions
F[0] = f1(y1, y2, y3, Y_old[0],Y_old[1],Y_old[2], dt)
F[1] = f2(y1, y2, y3, Y_old[0],Y_old[1],Y_old[2], dt)
F[2] = f3(y1, y2, y3, Y_old[0],Y_old[1],Y_old[2], dt)

# Computing New Values
Y_new = Y_old - alpha*(np.matmul(inv(J), F))

# Computing Max absolute error
error = np.max(np.abs(Y_new - Y_old))

# Updating old values
Y_old = 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
return [Y_new[0], Y_new[1], Y_new[2]]

def implicit_euler(y1, y2, y3, t_start, t_end, dt):
# Newton Rhapson Outer Integration Loop
t = np.arange(t_start,t_end,dt)
time = len(t)
Y_1 = np.zeros(time)
Y_2 = np.zeros(time)
Y_3 = np.zeros(time)

Y_1[0] = y1
Y_2[0] = y2
Y_3[0] = y3

y1_guess = 0
y2_guess = 0
y3_guess = 0

for i in range(1,time):
Y_1[i], Y_2[i], Y_3[i] = Newton_Rhapson(Y_1[i-1], Y_2[i-1], Y_3[i-1],y1_guess,y2_guess,y3_guess,dt)

# Final Results
y1_guess = Y_1[i]
y2_guess = Y_2[i]
y3_guess = Y_3[i]
return [t, Y_1, Y_2, Y_3]

t, y1, y2, y3 = implicit_euler(1,0,0,0,600,0.6)
plt.plot(t,y1)
plt.plot(t,y2)
plt.plot(t,y3)
plt.xlabel('Simulation Time')
plt.legend(['y1 graph','y2 graph','y3 graph'])
plt.show()

Observations:
In the above image, we see the approximation of the unknowns using the Multivariate Newton Rhapson iterative method. On the X-axis, we have the simulation time in seconds which is about 600 seconds (10 minutes). With respect to the given time, we see the change in values of the unknowns. As the approximations are implemented, the values change rapidly to different values. This behavior shows the Stiff characteristic of the given system. The step size was chosen for this simulation is 0.6. Even if the step size is changed, the system remains stable that is because we are using an Implicit method. Explicit methods are very sensitive to the time steps given which results in instability while implicit is less sensitive to the time steps.

Why?
The Stability Issue:
The principal reason for using implicit solution methods, which are more complex to program and require more computational effort in each solution step, is to allow for large time-step sizes. A simple qualitative model will help to illustrate how this works. Let x be a quantity whose value `x_(n+1)` we want to compute at time `t=(n+1)dt`, in terms of its value at time `t=n.dt`, i.e., `x_(n+1)=x_n + dt * S` where S is the rate of change in x.

In an explicit numerical method, S would be evaluated in terms of known quantities at the previous time step n.

As a general rule, it can be shown that when `alpha <=1` is very nearly equivalent to the stability condition for an explicit approximation. Another general rule is that the time-step sizes for explicit stability and accuracy are usually equivalent. Thus, when `alpha >=1`, an explicit method would be unstable, but implicit methods simply under-relax more to maintain the stability of the iterative solution. It is this increased damping, with the increase in time-step size, which produces inaccuracies in transient behavior.

Inference:
When the step size is sufficiently small only one iteration is necessary for convergence. This shows that the implicit formulation adds a smaller change to in one-time step than it would occur in an explicit method because of the under-relaxation factor that multiplies the time step. Explicit methods require a time-step size that limits the advance of the function approximation to less than one value per time step.