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.
void	exact_solution (const double &x, std::valarray< double > *y)
	Exact solution of the problem. Used for solution comparison.
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.
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.
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.

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]	x	independent variable
[in,out]	y	dependent variable

                                                             {
    y[0][0] = std::cos(x);
    y[0][1] = -std::sin(x);
}

Here is the call graph for this function:

main()

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

Main Function

                                 {
    double X0 = 0.f;                       /* initial value of x0 */
    double X_MAX = 10.F;                   /* upper limit of integration */
    std::valarray<double> Y0 = {1.f, 0.f}; /* initial value Y = y(x = x_0) */
    double step_size;
 
    if (argc == 1) {
        std::cout << "\nEnter the step size: ";
        std::cin >> step_size;
    } else {
        // use commandline argument as independent variable step size
        step_size = std::atof(argv[1]);
    }
 
    // get approximate solution
    double total_time = semi_implicit_euler(step_size, X0, X_MAX, &Y0, true);
    std::cout << "\tTime = " << total_time << " ms\n";
 
    /* compute exact solution for comparion */
    save_exact_solution(X0, X_MAX, step_size, Y0);
 
    return 0;
}

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]	x	independent variable(s)
[in,out]	y	dependent variable(s)
[in,out]	dy	first-derivative of dependent variable(s)

                                      {
    const double omega = 1.F;             // some const for the problem
    dy[0][0] = y[0][1];                   // x dot
    dy[0][1] = -omega * omega * y[0][0];  // y dot
}

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]	X0	initial value of independent variable
[in]	X_MAX	final value of independent variable
[in]	step_size	independent variable step size
[in]	Y0	initial values of dependent variables

                                                        {
    double x = X0;
    std::valarray<double> y = Y0;
 
    std::ofstream fp("exact.csv", std::ostream::out);
    if (!fp.is_open()) {
        std::perror("Error! ");
        return;
    }
    std::cout << "Finding exact solution\n";
 
    std::clock_t t1 = std::clock();
    do {
        fp << x << ",";
        for (int i = 0; i < y.size() - 1; i++) {
            fp << y[i] << ",";
        }
        fp << y[y.size() - 1] << "\n";
 
        exact_solution(x, &y);
 
        x += step_size;
    } while (x <= X_MAX);
 
    std::clock_t t2 = std::clock();
    double total_time = static_cast<double>(t2 - t1) / CLOCKS_PER_SEC;
    std::cout << "\tTime = " << total_time << " ms\n";
 
    fp.close();
}

Here is the call graph for this function: