From 990583f2c4e2aa298de9f77e28fb5613ffc9404b Mon Sep 17 00:00:00 2001
From: Krishna Vedala <7001608+kvedala@users.noreply.github.com>
Date: Wed, 10 Jun 2020 13:25:41 -0400
Subject: [PATCH] added semi-implicit euler method
---
numerical_methods/ode_semi_implicit_euler.c | 182 ++++++++++++++++++++
1 file changed, 182 insertions(+)
create mode 100644 numerical_methods/ode_semi_implicit_euler.c
diff --git a/numerical_methods/ode_semi_implicit_euler.c b/numerical_methods/ode_semi_implicit_euler.c
new file mode 100644
index 00000000..3375259b
--- /dev/null
+++ b/numerical_methods/ode_semi_implicit_euler.c
@@ -0,0 +1,182 @@
+/**
+ * \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
+ * [semi implicit Euler
+ * method](https://en.wikipedia.org/wiki/Semi-implicit_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 `semi_implicit_euler.csv`
+ * and the exact soltuion results in `exact.csv` for comparison. 
+ * \see ode_forward_euler.c, ode_midpoint_euler.c
+ */
+
+#include
+#include
+#include
+#include
+
+#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 midpoint-Euler
+ * method.
+ * @f[y_{n+1} = y_n + dx\, f\left(x_n+\frac{1}{2}dx,
+ * y_n + \frac{1}{2}dx\,f\left(x_n,y_n\right)\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$y_n+\frac{1}{2}dx\,f\left(x_n,y_n\right)@f$
+ */
+void semi_implicit_euler_step(double dx, double *x, double *y, double *dy)
+{
+ problem(x, y, dy);
+ double tmp_x = (*x) + 0.5 * dx;
+ double tmp_y[order];
+ int o;
+ for (o = 0; o < order; o++)
+ tmp_y[o] = y[o] + 0.5 * dx * dy[o];
+
+ problem(&tmp_x, tmp_y, dy);
+
+ for (o = 0; o < order; o++)
+ y[o] += dx * dy[o];
+}
+
+/**
+ * @brief Compute approximation using the midpoint-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 semi_implicit_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("semi_implicit_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
+ semi_implicit_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 = semi_implicit_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;
+}