diff --git a/doc/sphinx/source/examples/Time-Integrators/backward_euler.md b/doc/sphinx/source/examples/Time-Integrators/backward_euler.md index 7e383451..31860f91 100644 --- a/doc/sphinx/source/examples/Time-Integrators/backward_euler.md +++ b/doc/sphinx/source/examples/Time-Integrators/backward_euler.md @@ -6,7 +6,7 @@ $$ \frac{dy}{dt} = \sin^2(t) \cdot y $$ -with initial condition $y(0) = 0$ over the time interval $[0,5]$. +with initial condition $y(0) = 2$ over the time interval $[0,5]$. #### Mathematical Background @@ -26,16 +26,18 @@ where: #### Implementation Details Note that $f(t_{i+1}, y_{i+1})$ is not known due to $y_{i+1}$ being on both the left and right hand side. Therefore, we can solve for $y_{i+1}$ by finding the root: -$$ -y_{i+1} - y_i + h \cdot f(t_{i+1}, y_{i+1}) = 0 -$$ +$$y_{i+1} - y_i - h \cdot f(t_{i+1}, y_{i+1}) = 0$$ In the C++ example, we used the fixed-point iteration method for rootfinding due to it's simpler nature. #### Results -![Backward Euler gnuplot](../../_images/backward_euler.png) +
+ +
+Backward Euler solutions y(t) per time-step t. +
--- This example is implemented in: diff --git a/doc/sphinx/source/examples/Time-Integrators/figures/backward_euler.png b/doc/sphinx/source/examples/Time-Integrators/figures/backward_euler.png new file mode 100644 index 00000000..4d27d524 Binary files /dev/null and b/doc/sphinx/source/examples/Time-Integrators/figures/backward_euler.png differ diff --git a/doc/sphinx/source/examples/Time-Integrators/index.md b/doc/sphinx/source/examples/Time-Integrators/index.md index ea8c84a4..20223757 100644 --- a/doc/sphinx/source/examples/Time-Integrators/index.md +++ b/doc/sphinx/source/examples/Time-Integrators/index.md @@ -6,6 +6,7 @@ Time integration methods are essential for solving time-dependent PDEs, converti :maxdepth: 1 :caption: Contents +backward_euler RK2 RK4 ``` diff --git a/examples/cpp/backward_euler.cpp b/examples/cpp/backward_euler.cpp index 4dee7c96..16270533 100644 --- a/examples/cpp/backward_euler.cpp +++ b/examples/cpp/backward_euler.cpp @@ -1,9 +1,13 @@ /** * @file backward_euler.cpp - * @brief Solves a first-order ODE with backward Euler with an initial value condition. - * Equation: $\frac{dy}{dt} = \sin^2(t) \cdot y$ - * Domain: N/A - * Boundary Conditions: N/A + * @brief Solves a first-order ODE with backward Euler with an initial value + * condition. + * + * Equation: dy/dt = sin^2(t) * y + * + * Domain: N/A Boundary + * + * Conditions: N/A */ #include "mole.h" @@ -13,54 +17,63 @@ double f(double t, double y); int main() { - // Problem parameters - constexpr double h = 0.01; // Step-size h - arma::vec t = arma::regspace(0, h, 5); // Calculates up to y(5) at step-size h - arma::mat y = arma::mat(1, t.size(), arma::fill::zeros); - constexpr double tol = 1e-6; // tolerance for fixed-point iteration + // Problem parameters + constexpr double h = 0.01; // Step-size h + arma::vec t = arma::regspace(0, h, 5); // Calculates up to y(5) at step-size h + arma::mat y = arma::mat(1, t.size(), arma::fill::zeros); + constexpr double tol = 1e-6; // Tolerance for fixed-point iteration + constexpr int max_iter = 100; // Max iterations for fixed-point iteration - // Initial conditions - y(0) = 2.0; - - // Backward Euler - for (int i = 0; i < t.size()-1; i++) { - double y_old = y(i); - for (int k = 0; k < 100; k++) { // fixed-point iteration for rootfinding - y(i+1) = y_old + h*f(t(i+1), y_old); // Backward Euler - if (std::abs(y(i+1) - y_old) < tol) { // Stopping criteria for approximate relative error - break; - } - } - } + // Initial conditions + y(0) = 2.0; - // Create a GNUplot script file - std::ofstream plot_script("plot.gnu"); - if (!plot_script) { - std::cerr << "Error: Failed to create GNUplot script.\n"; - return 1; + // Backward Euler + for (int i = 0; i < t.size() - 1; i++) { + double y_old = y(i); + double xn = y_old; // Initial input for fixed-point iteration + double xnp1; // x at iteration n+1 + for (int n = 0; n < max_iter; + n++) { // fixed-point iteration for rootfinding + xnp1 = y_old + h * f(t(i + 1), xn); // Backward Euler + if (std::abs(xnp1 - xn) < + tol) { // Stopping criteria for approximate relative error + break; + } + xn = xnp1; } - plot_script << "set title 'Approximation to y(t) using Backward Euler'\n"; - plot_script << "set xlabel 't'\n"; - plot_script << "set ylabel 'y'\n"; - plot_script << "plot '-' using 1:2 with lines\n"; - - // Print the time and solution values - for (int i = 0; i < t.size(); ++i) { - // output to stdout - std::cout << t(i) << " " << y(i) << "\n"; - // AND output to plot_script (plot.gnu) - plot_script << t(i) << " " << y(i) << "\n"; + y(i + 1) = xnp1; // root found + if (std::abs(xnp1 - xn) > tol) { // Fixed-point did not converge + std::cerr << "Warning: Fixed-point iteration did not converge. Try " + "reducing the step-size or adjusting the tolerance."; } - plot_script.close(); + } - // Execute GNUplot using the script - if (system("gnuplot -persist plot.gnu") != 0) { - std::cerr << "Error: Failed to execute GNUplot.\n"; - return 1; - } + // Create a GNUplot script file + std::ofstream plot_script("plot.gnu"); + if (!plot_script) { + std::cerr << "Error: Failed to create GNUplot script.\n"; + return 1; + } + plot_script << "set title 'Approximation to y(t) using Backward Euler'\n"; + plot_script << "set xlabel 't'\n"; + plot_script << "set ylabel 'y'\n"; + plot_script << "plot '-' using 1:2 with lines\n"; + + // Print the time and solution values + for (int i = 0; i < t.size(); ++i) { + // output to stdout + std::cout << t(i) << " " << y(i) << "\n"; + // AND output to plot_script (plot.gnu) + plot_script << t(i) << " " << y(i) << "\n"; + } + plot_script.close(); + + // Execute GNUplot using the script + if (system("gnuplot -persist plot.gnu") != 0) { + std::cerr << "Error: Failed to execute GNUplot.\n"; + return 1; + } } // Function definition for f(t, y) -double f(double t, double y) { - return std::pow(std::sin(t), 2) * y; -} \ No newline at end of file +double f(double t, double y) { return std::pow(std::sin(t), 2) * y; } \ No newline at end of file