Multivariate Newton Rhapson Solver

OBJECTIVE:-

In the current challenge, we are going to be programming a multivariate Newton raphson solver for the following system of equation

In the given equation system we can say they are combined and non-linear, and we can verify their association y1,y2,y3 can only be solved the three equations together. They say it\'s non-linear because in the equations the unknown variable is multiplied by itself or by each other.

Here, for each time phase set for the given time of simulation, the Newton raphson method is used. The difference is that, as there is more than one unknown, the matrix approximation in multivariate form is determined.

The general principle of the Multivariate newton raphson is given below

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 which contains the given system of equation

The matrix representation of the above principle is given below

Alpha is multiplied in the above figure to the jacobian reverse and the word matrix of functions.

We determined the Jacobian matrix numerically for the present task since analytical derivatives aren\'t always available.

The principle of the numerical Jacobian is expressed below.

In this scenario, the simulation is running for 10 minutes. The relaxing factor is taken as 1 in the newton raphson solver.

We are also using Implicit Euler or otherwise know as a backward differencing method for the outer time integration loop for this challenge.

The backward Differencing method is a finite approximation method that uses the function values at (x) and (x-h) where h is the step size. The general principle for the method is mentioned below.

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 condition for the unknowns are

y₁ = 1,

y₂=0,

y₃=0.

Here we are initializing the ODE’s in the form

The general layout of the code is given below

• Outer time integration loop
• Assembling the ODE system
• Newton raphson solver till convergence
• Get new values and update them
• End of the outer loop.

The code for the solver is given below:

\"\"\"
Pogram to solve system of coupled non-linear equations using multivariate Newton Rhapson 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

#intializtion
t_initial = 0
#t_final = 10*60
t_final = 600
dt = 0.05


#intial 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):

		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))

		error = 9e9
		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 function
				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 the new values
				Y_new = Y_old - alpha*(np.matmul(inv(J), F))

				#Computing the 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_initial, t_final, dt):

				t = np.arange(t_initial,t_final,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 result

						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.05)

plt.plot(t,y1,color = \'red\')
plt.plot(t,y2,color = \'green\')
plt.plot(t,y3,color = \'blue\')
plt.grid(\'on\')
plt.xlabel(\'Simulation time\')
plt.ylabel(\'Initial guess values\')
plt.legend([\'y1\',\'y2\',\'y3\'])
plt.savefig(\'Multi-varable-newton-rhapson\')
plt.show()

Result:

#intializtion
t_initial = 0
#t_final = 10*60
t_final = 600
dt = 0.5

Result:

Observation:

As the solution marches forward in time, the values change rapidly to different values. This behavior shows the stiff characteristic of the system.

As for the value of the time step,0.05 is used for solving this problem but changing the time step does not seem to affect the results. thought if we give a small value of time step the solution marches more forward in time as the end time remains the same and for a large value of time step, the solution still remains stable and is able to converge and provide meaning results.

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.

Though explicit schemes are easy to solve compared to the implicit schemes.