Multivariate Newton Rhapson Solver

In this project, coupled Ordinary Differential Equations are solved by using Multivariate Newton Rhapson solver. The equations that are to solved is represented below:

`(dy_1)/dt = (-0.04xxy_1) + (10^4xxy_2xxy_3)`

`(dy_2)/dt = (0.04xxy_1) - (10^4xxy_2xxy_3) - (3xx10^7xxy_2^2)`

`(dy_3)/dt = 3xx10^7xxy_2^2`

Now, the ODEs are converted into the following functions to find solution, these functions are mentioned below:

`y_1^(n-1) + dtxx(-0.04xxy_1^n + 10^4xxy_2^nxxy_3^n) - y_1^n = 0`

`y_2^(n-1) + dtxx(0.04xxy_1^n - 10^4xxy_2^nxxy_3^n - 3xx10^7xxy_2^nxxy_2^n)-y_2^n = 0`

`y_3^(n-1) + dtxx(3xx10^7xxy_2^nxxy_2^n) - y_3^n = 0`

In the above equations usied for computing function, `y_x^(n-1)` denotes old time value or initial value in case of computing the first time step value whereas `y_x^n` denotes the new time value that is to be computed. The `dt`here is the time step size. Since only old time values or initial value is used in computing all other new time step values through iteration the method used here is the implicit method also know as the backward difference method.

Jacobi matrix is a matrix containing all types of possible derivates of the three function mentioned above. Here the Jacobi matrix is found by using numerical method using functions with step size h of 1e-8. The derivates of the functions are computated numerically as follows:

`dy_x/dt = (f(y_x + h) - f(y_x))/h`

Initial values are taken as follow:

y1 = 1

y2 = 0

y3 = 0

Relaxation factor `alpha` is taken as 1 for the solver and a time step size of 1 is chosen with tolerance kept as 1e-6 in the program and the initial error set as 9e9. The computation is carried out for 10 minutes i.e., 600 seconds.

The formulae used for calculating the new iterative value at that particular time is given below:

`y_(\"new iter\") = y_(\"old iter\") - alphaxx(J^-1xxF)`

here,

`y_(\"new iter\")` is the new value of iteration at that particular time

`y_(\"old iter\")` is the old value of iteration at the particular time

`alpha` is the relaxation factor

`J^(-1)` is the inverse Jacobi matrix

`F` is the column matrix consisting of the three functions made up of the equation shown before.

The code of the program is attached below:

https://drive.google.com/file/d/1hXQGCgf0WZuNp2RdzGMJ2JVodHhVm6yd/view?usp=sharing

The results obtained are mapped in graph plotting values of the dependent variables `y_1`,`y_2` and `y_3` against time. The image of the graph obtained for dt = 0.2s is pasted below:

Since, the method used is essentially implicit, which is unconditionally stable for various time step sizes the solution obtained is stable but not accurate. The only way to get accurate solution is by reducing the time step size till the values do not change significantly. By following this method the ideal time step size is found to be 0.2s. The results from which the conclusion was made are given below.

For dt = 0.25s, the results were:

For dt = 0.2s, the results were:

For dt = 0.15s, the results were:

Finally,the results for dt = 0.1s were: