Week 6 - Multivariate Newton Rhapson Solver

                                                   MULTIVARIATE NEWTON RAPHSON SOLVER FOR ODE'S

Objective

Solve the problem using Implicit Euler Method/Backward Differencing, assume the initial conditions as 1 for y1 and 0 for y2 and y3 respectively.

Theory

http://fourier.eng.hmc.edu/e176/lectures/NM/node21.html

The equations are independent and coupled in nature.we have 3 non linearcoupled equation with 3 unknowns inorder to solve this we use newton raphson but since we have multiple equation we use multivariate newton raphson scheme.

`x^(k+1)=x^(k)-alpha*J^(-1)*F`

where`x^k` is old value

`x^(k+1)` new value guessing from old value

alpha is the relaxation factor

J is the jacobian matrix

F is the function " contains system of equation"

Code:

import numpy as np
from numpy.linalg import inv
import matplotlib.pyplot as plt

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

def f2(y2_old,y1,y2,y3,dt):
    return(y2_old)+(dt*((0.04*y1)-(10**4*y2*y3)-(3*10**7*y2**2)))-y2

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

#defining jacobian matrix  (backward defferencing)

def jacobian (y1_old,y2_old,y3_old,y1,y2,y3,dt):
    J= np.ones ((3,3))

#specifing the time step for solving the newtonraphson method

    h = 1e-6

#row1

    J[0,0] = (f1(y1_old,y1+h,y2,y3,dt)-f1(y1_old,y1,y2,y3,dt))/h

    J[0,1] = (f1(y1_old,y1,y2+h,y3,dt)-f1(y1_old,y1,y2,y3,dt))/h

    J[0,2] = (f1(y1_old,y1,y2,y3+h,dt)-f1(y1_old,y1,y2,y3,dt))/h
    

#row2
    J[1,0]=(f2(y2_old,y1+h,y2,y3,dt)-f3(y2_old,y1,y2,y3,dt))/h
    J[1,1] = (f2(y2_old,y1,y2+h,y3,dt)-f3(y2_old,y1,y2,y3,dt))/h 
    J[1,2] =(f2(y2_old,y1,y2,y3+h,dt)-f3(y2_old,y1,y2,y3,dt))/h
    
#row 3    
    
    J[2,0]=(f3(y3_old,y1+h,y2,y3,dt)-f3(y3_old,y1,y2,y3,dt))/h
    J[2,1] = (f3(y3_old,y1,y2+h,y3,dt)-f3(y3_old,y1,y2,y3,dt))/h 
    J[2,2] =(f3(y3_old,y1,y2,y3+h,dt)-f3(y3_old,y1,y2,y3,dt))/h



    return J

# storing old values

y_old = np.ones ((3,1))
y_new = np.ones ((3,1))
y_guess = np.ones ((3,1))

y1 = 1
y2 = 0
y3 = 0

#specifying the time step

dt =0.05
# running for 600seconds with dt value 0.05
t = np.arange(0, 600, dt)
Num_of_iter = len(t)
y = np.ones ((3, Num_of_iter))


F = np.copy(y_guess)
alpha = 1
tol = 1e-8

iter = 0
a=0

for i in range (0, Num_of_iter):
    y_old[0] = y1

    y_old [1] = y2

    y_old [2] = y3

    y_guess [0] = y1

    y_guess [1] = y2

    y_guess [2] = y3
    error = 9e9


    
    while error>tol:
        
        
        

# updating the F vector for every time step

        F[0] = f1(y_old [0],y_guess [0],y_guess [1],y_guess [2],dt);

        F[1] = f2(y_old[1],y_guess [0],y_guess [1],y_guess [2], dt);

        F[2] = f3(y_old[2],y_guess [0],y_guess [1],y_guess [2], dt);

        J = jacobian (y_old[0],y_old[1],y_old [2],y_guess [0],y_guess[1],y_guess [2], dt)

        y_new = y_guess-alpha*np.matmul(inv(J), F) 
        error = np.max(np.abs(y_new-y_guess))

        y_guess = y_new

        log_message = 'iteration# = {0} y1={1} y2={2} y3={3} error ={4}'. format (iter,y_new[0],y_new[1],y_new[2], error)

        print (log_message)

        iter = iter+1

    y1 = y_new[0]

    y2 = y_new[1]

    y3 = y_new[2]
   
    y[0,a] =y_new[0]
    y[1,a] =y_new[1]
    y[2,a] = y_new[2]

    a=a+1

plt.plot(y[0])

plt.plot(y[1])

plt.plot(y[2])
plt.xlabel('Number of steps')

plt.ylabel('Variables (Species)')

plt.grid ('on')

plt.title('Multivariate Newton Raphson Method for dt = 0.05') 

plt.legend (['y1', 'y2', 'y3'])

plt.show()

Observation:

for dt=0.05

for dt=0.5

Comparing both the plots first observation we find is that ,behaviour of equation y2 different from y1 and y2.Although we can't observe the change,value of y2 is really small which is not shown in the scale of the graph.The values changes rapidly with respect to time shows the behaviour of Stiff Characteristics of the system.

Now comparing the time steps,changing the time steps doesn't appear to change the results,by using smaller time step we get increased iteration and computational time.In the end the solution still remains stable and was able to converge.This is because we have used an implicit scheme ,which is stable irrespective of timesteps.Although generaly since explicit schemes are dependent on timesteps,comparetively explicit schemes are easy to solve than implicit.