## C, , ode_forward_euler.c

/**
* \file
* \authors [Krishna Vedala](https://github.com/kvedala)
* \brief Solve a multivariable first order [ordinary differential equation
* (ODEs)](https://en.wikipedia.org/wiki/Ordinary_differential_equation) using
* [forward Euler
* method](https://en.wikipedia.org/wiki/Numerical_methods_for_ordinary_differential_equations#Euler_method)
*
* \details
* The ODE being solved is:
* \f{eqnarray*}{
* \dot{u} &=& v\\
* \dot{v} &=& -\omega^2 u\\
* \omega &=& 1\\
* [x_0, u_0, v_0] &=& [0,1,0]\qquad\ldots\text{(initial values)}
* \f}
* The exact solution for the above problem is:
* \f{eqnarray*}{
* u(x) &=& \cos(x)\\
* v(x) &=& -\sin(x)\\
* \f}
* The computation results are stored to a text file forward_euler.csv and the
* exact soltuion results in exact.csv for comparison.
* <img
* src="https://raw.githubusercontent.com/TheAlgorithms/C/docs/images/numerical_methods/ode_forward_euler.svg"
* alt="Implementation solution"/>
*
* To implement [Van der Pol
* oscillator](https://en.wikipedia.org/wiki/Van_der_Pol_oscillator), change the
* ::problem function to:
* cpp
* const double mu = 2.0;
* dy[0] = y[1];
* dy[1] = mu * (1.f - y[0] * y[0]) * y[1] - y[0];
* 
* \see ode_midpoint_euler.c, ode_semi_implicit_euler.c
*/

#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <time.h>

#define order 2 /**< number of dependent variables in ::problem */

/**
* @brief Problem statement for a system with first-order differential
* equations. Updates the system differential variables.
* \note This function can be updated to and ode of any order.
*
* @param[in] 		x 		independent variable(s)
* @param[in,out]	y		dependent variable(s)
* @param[in,out]	dy	    first-derivative of dependent variable(s)
*/
void problem(const double *x, double *y, double *dy)
{
const double omega = 1.F;       // some const for the problem
dy[0] = y[1];                   // x dot
dy[1] = -omega * omega * y[0];  // y dot
}

/**
* @brief Exact solution of the problem. Used for solution comparison.
*
* @param[in] 		x 		independent variable
* @param[in,out]	y		dependent variable
*/
void exact_solution(const double *x, double *y)
{
y[0] = cos(x[0]);
y[1] = -sin(x[0]);
}

/**
* @brief Compute next step approximation using the forward-Euler
* method. @f[y_{n+1}=y_n + dx\cdot f\left(x_n,y_n\right)@f]
* @param[in] 		dx	step size
* @param[in,out] 	x	take \f$x_n\f$ and compute \f$x_{n+1}\f$
* @param[in,out] 	y	take \f$y_n\f$ and compute \f$y_{n+1}\f$
* @param[in,out]	dy	compute \f$f\left(x_n,y_n\right)\f$
*/
void forward_euler_step(const double dx, const double *x, double *y, double *dy)
{
int o;
problem(x, y, dy);
for (o = 0; o < order; o++) y[o] += dx * dy[o];
}

/**
* @brief Compute approximation using the forward-Euler
* method in the given limits.
* @param[in] 		dx  	step size
* @param[in]   	x0  	initial value of independent variable
* @param[in] 	    x_max	final value of independent variable
* @param[in,out] 	y	    take \f$y_n\f$ and compute \f$y_{n+1}\f$
* @param[in] save_to_file	flag to save results to a CSV file (1) or not (0)
* @returns time taken for computation in seconds
*/
double forward_euler(double dx, double x0, double x_max, double *y,
char save_to_file)
{
double dy[order];

FILE *fp = NULL;
if (save_to_file)
{
fp = fopen("forward_euler.csv", "w+");
if (fp == NULL)
{
perror("Error! ");
return -1;
}
}

/* start integration */
clock_t t1 = clock();
double x = x0;
do  // iterate for each step of independent variable
{
if (save_to_file && fp)
fprintf(fp, "%.4g,%.4g,%.4g\n", x, y[0], y[1]);  // write to file
forward_euler_step(dx, &x, y, dy);  // perform integration
x += dx;                            // update step
} while (x <= x_max);  // till upper limit of independent variable
/* end of integration */
clock_t t2 = clock();

if (save_to_file && fp)
fclose(fp);

return (double)(t2 - t1) / CLOCKS_PER_SEC;
}

/**
Main Function
*/
int main(int argc, char *argv[])
{
double X0 = 0.f;          /* initial value of x0 */
double X_MAX = 10.F;      /* upper limit of integration */
double Y0[] = {1.f, 0.f}; /* initial value Y = y(x = x_0) */
double step_size;

if (argc == 1)
{
printf("\nEnter the step size: ");
scanf("%lg", &step_size);
}
else
// use commandline argument as independent variable step size
step_size = atof(argv[1]);

// get approximate solution
double total_time = forward_euler(step_size, X0, X_MAX, Y0, 1);
printf("\tTime = %.6g ms\n", total_time);

/* compute exact solution for comparion */
FILE *fp = fopen("exact.csv", "w+");
if (fp == NULL)
{
perror("Error! ");
return -1;
}
double x = X0;
double *y = &(Y0[0]);
printf("Finding exact solution\n");
clock_t t1 = clock();

do
{
fprintf(fp, "%.4g,%.4g,%.4g\n", x, y[0], y[1]);  // write to file
exact_solution(&x, y);
x += step_size;
} while (x <= X_MAX);

clock_t t2 = clock();
total_time = (t2 - t1) / CLOCKS_PER_SEC;
printf("\tTime = %.6g ms\n", total_time);
fclose(fp);

return 0;
}