Algorithms_in_C++  1.0.0
Set of algorithms implemented in C++.
ode_semi_implicit_euler.cpp File Reference

Solve a multivariable first order ordinary differential equation (ODEs) using semi implicit Euler method More...

#include <cmath>
#include <ctime>
#include <fstream>
#include <iostream>
#include <valarray>
Include dependency graph for ode_semi_implicit_euler.cpp:

Functions

void problem (const double &x, std::valarray< double > *y, std::valarray< double > *dy)
 Problem statement for a system with first-order differential equations. Updates the system differential variables. More...
 
void exact_solution (const double &x, std::valarray< double > *y)
 Exact solution of the problem. Used for solution comparison. More...
 
void semi_implicit_euler_step (const double dx, const double &x, std::valarray< double > *y, std::valarray< double > *dy)
 Compute next step approximation using the semi-implicit-Euler method. More...
 
double semi_implicit_euler (double dx, double x0, double x_max, std::valarray< double > *y, bool save_to_file=false)
 Compute approximation using the semi-implicit-Euler method in the given limits. More...
 
void save_exact_solution (const double &X0, const double &X_MAX, const double &step_size, const std::valarray< double > &Y0)
 
int main (int argc, char *argv[])
 

Detailed Description

Solve a multivariable first order ordinary differential equation (ODEs) using semi implicit Euler method

Authors
Krishna Vedala

The ODE being solved is:

\begin{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)} \end{eqnarray*}

The exact solution for the above problem is:

\begin{eqnarray*} u(x) &=& \cos(x)\\ v(x) &=& -\sin(x)\\ \end{eqnarray*}

The computation results are stored to a text file semi_implicit_euler.csv and the exact soltuion results in exact.csv for comparison. Implementation solution

To implement Van der Pol oscillator, change the problem function to:

const double mu = 2.0;
dy[0] = y[1];
dy[1] = mu * (1.f - y[0] * y[0]) * y[1] - y[0];
See also
ode_midpoint_euler.cpp, ode_forward_euler.cpp

Function Documentation

◆ exact_solution()

void exact_solution ( const double &  x,
std::valarray< double > *  y 
)

Exact solution of the problem. Used for solution comparison.

Parameters
[in]xindependent variable
[in,out]ydependent variable
66  {
67  y[0][0] = std::cos(x);
68  y[0][1] = -std::sin(x);
69 }
Here is the call graph for this function:

◆ main()

int main ( int  argc,
char *  argv[] 
)

Main Function

189  {
190  double X0 = 0.f; /* initial value of x0 */
191  double X_MAX = 10.F; /* upper limit of integration */
192  std::valarray<double> Y0 = {1.f, 0.f}; /* initial value Y = y(x = x_0) */
193  double step_size;
194 
195  if (argc == 1) {
196  std::cout << "\nEnter the step size: ";
197  std::cin >> step_size;
198  } else {
199  // use commandline argument as independent variable step size
200  step_size = std::atof(argv[1]);
201  }
202 
203  // get approximate solution
204  double total_time = semi_implicit_euler(step_size, X0, X_MAX, &Y0, true);
205  std::cout << "\tTime = " << total_time << " ms\n";
206 
207  /* compute exact solution for comparion */
208  save_exact_solution(X0, X_MAX, step_size, Y0);
209 
210  return 0;
211 }
Here is the call graph for this function:

◆ problem()

void problem ( const double &  x,
std::valarray< double > *  y,
std::valarray< double > *  dy 
)

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.
Parameters
[in]xindependent variable(s)
[in,out]ydependent variable(s)
[in,out]dyfirst-derivative of dependent variable(s)
54  {
55  const double omega = 1.F; // some const for the problem
56  dy[0][0] = y[0][1]; // x dot
57  dy[0][1] = -omega * omega * y[0][0]; // y dot
58 }

◆ save_exact_solution()

void save_exact_solution ( const double &  X0,
const double &  X_MAX,
const double &  step_size,
const std::valarray< double > &  Y0 
)

Function to compute and save exact solution for comparison

Parameters
[in]X0initial value of independent variable
[in]X_MAXfinal value of independent variable
[in]step_sizeindependent variable step size
[in]Y0initial values of dependent variables
155  {
156  double x = X0;
157  std::valarray<double> y = Y0;
158 
159  std::ofstream fp("exact.csv", std::ostream::out);
160  if (!fp.is_open()) {
161  std::perror("Error! ");
162  return;
163  }
164  std::cout << "Finding exact solution\n";
165 
166  std::clock_t t1 = std::clock();
167  do {
168  fp << x << ",";
169  for (int i = 0; i < y.size() - 1; i++) {
170  fp << y[i] << ",";
171  }
172  fp << y[y.size() - 1] << "\n";
173 
174  exact_solution(x, &y);
175 
176  x += step_size;
177  } while (x <= X_MAX);
178 
179  std::clock_t t2 = std::clock();
180  double total_time = static_cast<double>(t2 - t1) / CLOCKS_PER_SEC;
181  std::cout << "\tTime = " << total_time << " ms\n";
182 
183  fp.close();
184 }
Here is the call graph for this function:
semi_implicit_euler
double semi_implicit_euler(double dx, double x0, double x_max, std::valarray< double > *y, bool save_to_file=false)
Compute approximation using the semi-implicit-Euler method in the given limits.
Definition: ode_semi_implicit_euler.cpp:103
std::clock_t
std::cos
T cos(T... args)
std::atof
T atof(T... args)
std::clock
T clock(T... args)
std::cout
std::ofstream
STL class.
std::perror
T perror(T... args)
std::valarray< double >
std::sin
T sin(T... args)
save_exact_solution
void save_exact_solution(const double &X0, const double &X_MAX, const double &step_size, const std::valarray< double > &Y0)
Definition: ode_semi_implicit_euler.cpp:153
exact_solution
void exact_solution(const double &x, std::valarray< double > *y)
Exact solution of the problem. Used for solution comparison.
Definition: ode_semi_implicit_euler.cpp:66
std::cin