mirror of https://github.com/TheAlgorithms/C
Merge pull request #15 from kvedala/forward-euler
ODE solver using euler methods
This commit is contained in:
Krishna Vedala
GitHub
commit
8fccc20524
|
|
@ -15,7 +15,6 @@ Contributors guide: https://github.com/TheAlgorithms/C-Plus-Plus/CONTRIBUTING.md
|
|||
- [ ] Relevant documentation/comments is changed or added
|
||||
- [ ] PR title follows semantic [commit guidelines](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/CONTRIBUTING.md#Commit-Guidelines)
|
||||
- [ ] Search previous suggestions before making a new one, as yours may be a duplicate.
|
||||
- [ ] Sort by alphabetical order
|
||||
- [ ] I acknowledge that all my contributions will be made under the project's license.
|
||||
|
||||
Notes: <!-- Please add a one-line description for developers or pull request viewers -->
|
||||
|
|
@ -246,6 +246,9 @@
|
|||
* [Mean](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/numerical_methods/mean.c)
|
||||
* [Median](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/numerical_methods/median.c)
|
||||
* [Newton Raphson Root](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/numerical_methods/newton_raphson_root.c)
|
||||
* [Ode Forward Euler](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/numerical_methods/ode_forward_euler.c)
|
||||
* [Ode Midpoint Euler](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/numerical_methods/ode_midpoint_euler.c)
|
||||
* [Ode Semi Implicit Euler](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/numerical_methods/ode_semi_implicit_euler.c)
|
||||
* [Qr Decompose](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/numerical_methods/qr_decompose.h)
|
||||
* [Qr Decomposition](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/numerical_methods/qr_decomposition.c)
|
||||
* [Qr Eigen Values](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/numerical_methods/qr_eigen_values.c)
|
||||
|
|
|
|||
|
|
@ -24,11 +24,14 @@
|
|||
* \n Steps:
|
||||
* 1. `r1 = rand() % 100` gets a random number between 0 and 99
|
||||
* 2. `r2 = r1 / 100` converts random number to be between 0 and 0.99
|
||||
* 3. scale and offset the random number to given range of \f$[a,b]\f$
|
||||
* 3. scale and offset the random number to given range of \f$[a,b)\f$
|
||||
* \f[
|
||||
* y = (b - a) \times \frac{\text{(random number between 0 and RAND_MAX)} \;
|
||||
* \text{mod}\; 100}{100} + a \f]
|
||||
*
|
||||
* \param[in] a lower limit
|
||||
* \param[in] b upper limit
|
||||
* \returns random number in the range \f$[a,b]\f$
|
||||
* \returns random number in the range \f$[a,b)\f$
|
||||
*/
|
||||
double _random(double a, double b)
|
||||
{
|
||||
|
|
@ -184,7 +187,12 @@ void kohonen_som_tracer(double **X, double *const *W, int num_samples,
|
|||
/** Creates a random set of points distributed *near* the circumference
|
||||
* of a circle and trains an SOM that finds that circular pattern. The
|
||||
* generating function is
|
||||
* \f{eqnarray*}{ \f}
|
||||
* \f{eqnarray*}{
|
||||
* r &\in& [1-\delta r, 1+\delta r)\\
|
||||
* \theta &\in& [0, 2\pi)\\
|
||||
* x &=& r\cos\theta\\
|
||||
* y &=& r\sin\theta
|
||||
* \f}
|
||||
*
|
||||
* \param[out] data matrix to store data in
|
||||
* \param[in] N number of points required
|
||||
|
|
@ -269,8 +277,16 @@ void test1()
|
|||
|
||||
/** Creates a random set of points distributed *near* the locus
|
||||
* of the [Lamniscate of
|
||||
* Gerono](https://en.wikipedia.org/wiki/Lemniscate_of_Gerono) and trains an SOM
|
||||
* that finds that circular pattern. \param[out] data matrix to store data in
|
||||
* Gerono](https://en.wikipedia.org/wiki/Lemniscate_of_Gerono).
|
||||
* \f{eqnarray*}{
|
||||
* \delta r &=& 0.2\\
|
||||
* \delta x &\in& [-\delta r, \delta r)\\
|
||||
* \delta y &\in& [-\delta r, \delta r)\\
|
||||
* \theta &\in& [0, \pi)\\
|
||||
* x &=& \delta x + \cos\theta\\
|
||||
* y &=& \delta y + \frac{\sin(2\theta)}{2}
|
||||
* \f}
|
||||
* \param[out] data matrix to store data in
|
||||
* \param[in] N number of points required
|
||||
*/
|
||||
void test_lamniscate(double *const *data, int N)
|
||||
|
|
@ -354,10 +370,14 @@ void test2()
|
|||
free(W);
|
||||
}
|
||||
|
||||
/** Creates a random set of points distributed *near* the locus
|
||||
* of the [Lamniscate of
|
||||
* Gerono](https://en.wikipedia.org/wiki/Lemniscate_of_Gerono) and trains an SOM
|
||||
* that finds that circular pattern. \param[out] data matrix to store data in
|
||||
/** Creates a random set of points distributed in four clusters in
|
||||
* 3D space with centroids at the points
|
||||
* * \f$(0,5, 0.5, 0.5)\f$
|
||||
* * \f$(0,5,-0.5, -0.5)\f$
|
||||
* * \f$(-0,5, 0.5, 0.5)\f$
|
||||
* * \f$(-0,5,-0.5, -0.5)\f$
|
||||
*
|
||||
* \param[out] data matrix to store data in
|
||||
* \param[in] N number of points required
|
||||
*/
|
||||
void test_3d_classes(double *const *data, int N)
|
||||
|
|
|
|||
|
|
@ -0,0 +1,183 @@
|
|||
/**
|
||||
* \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/kvedala/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;
|
||||
}
|
||||
|
|
@ -0,0 +1,191 @@
|
|||
/**
|
||||
* \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
|
||||
* [midpoint Euler
|
||||
* method](https://en.wikipedia.org/wiki/Midpoint_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 `midpoint_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_midpoint_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_forward_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 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 midpoint_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 midpoint_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("midpoint_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
|
||||
midpoint_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 = midpoint_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;
|
||||
}
|
||||
|
|
@ -0,0 +1,193 @@
|
|||
/**
|
||||
* \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)
|
||||
*
|
||||
* \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 `semi_implicit_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_semi_implicit_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];
|
||||
* ```
|
||||
* <a href="https://en.wikipedia.org/wiki/Van_der_Pol_oscillator"><img
|
||||
* src="https://raw.githubusercontent.com/kvedala/C/docs/images/numerical_methods/van_der_pol_implicit_euler.svg"
|
||||
* alt="Van der Pol Oscillator solution"/></a>
|
||||
*
|
||||
* \see ode_forward_euler.c, ode_midpoint_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 semi-implicit-Euler
|
||||
* method.
|
||||
* @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;
|
||||
}
|
||||
Loading…
Reference in New Issue