added first order ODE solver using Euler method

numerical_methods/ode_forward_euler.c

@@ -0,0 +1,151 @@
+/**
+ * \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)
+ *
+ * \description
+ * 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/kvedala/C/docs/images/numerical_methods/ode_forward_euler.svg"
+ * alt="Implementation solution"/>
+ */
+
+#include <math.h>
+#include <stdio.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(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(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$
+ * @param[in]    order  order of algorithm implementation
+ */
+void forward_euler(double dx, double *x, double *y, double *dy)
+{
+    int o;
+    problem(x, y, dy);
+    for (o = 0; o < order; o++)
+        y[o] += dx * dy[o];
+    *x += dx;
+}
+
+/**
+    Main Function
+*/
+int main(int argc, char *argv[])
+{
+    double X0 = 0.f;          /* initial value of f(x = x0) */
+    double Y0[] = {1.f, 0.f}; /* initial value Y = y(x = x_0) */
+
+    double dx, dy[order];
+    double x = X0, *y = &(Y0[0]);
+    double X_MAX = 10.F; /* upper limit of integration */
+
+    if (argc == 1)
+    {
+        printf("\nEnter the step size: ");
+        scanf("%lg", &dx);
+    }
+    else
+        // use commandline argument as independent variable step size
+        dx = atof(argv[1]);
+
+    clock_t t1, t2;
+    double total_time;
+
+    FILE *fp = fopen("forward_euler.csv", "w+");
+    if (fp == NULL)
+    {
+        perror("Error! ");
+        return -1;
+    }
+    printf("Computing using 'Forward Euler' algorithm\n");
+
+    /* start integration */
+    t1 = clock();
+    do // iterate for each step of independent variable
+    {
+        fprintf(fp, "%.4g,%.4g,%.4g\n", x, y[0], y[1]); // write to file
+        forward_euler(dx, &x, y, dy);                   // perform integration
+    } while (x <= X_MAX); // till upper limit of independent variable
+    /* end of integration */
+
+    t2 = clock();
+    fclose(fp);
+
+    total_time = (t2 - t1) / CLOCKS_PER_SEC;
+    printf("\tTime taken = %.6g ms\n", total_time);
+
+    /* compute exact solution for comparion */
+    fp = fopen("exact.csv", "w+");
+    if (fp == NULL)
+    {
+        perror("Error! ");
+        return -1;
+    }
+    x = X0;
+    y = Y0;
+    printf("Finding exact solution\n");
+    t1 = clock();
+    do
+    {
+        fprintf(fp, "%.4g,%.4g,%.4g\n", x, y[0], y[1]); // write to file
+        exact_solution(&x, y);
+        x += dx;
+    } while (x <= X_MAX);
+    t2 = clock();
+    total_time = (t2 - t1) / CLOCKS_PER_SEC;
+    printf("\tTime = %.6g ms\n", total_time);
+    fclose(fp);
+
+    return 0;
+}