Correction to Backward Euler C++ Example and Documentation #249
Conversation
jbrzensk
left a comment
jbrzensk
left a comment
There was a problem hiding this comment.
Good catch! missed the "fzero" in Matlab/Octave, thats the solve for that method.
The output also matches very well:
4.85 24.5469
4.86 24.7894
4.87 25.0336
4.88 25.2793
4.89 25.5266
4.9 25.7754
4.91 26.0256
4.92 26.2772
4.93 26.5302
4.94 26.7844
4.95 27.0398
4.96 27.2964
4.97 27.554
4.98 27.8127
4.99 28.0723
5 28.3329
versus from matlab:
Time: 4.850, Val: 24.547
Time: 4.860, Val: 24.789
Time: 4.870, Val: 25.034
Time: 4.880, Val: 25.279
Time: 4.890, Val: 25.527
Time: 4.900, Val: 25.775
Time: 4.910, Val: 26.026
Time: 4.920, Val: 26.277
Time: 4.930, Val: 26.530
Time: 4.940, Val: 26.784
Time: 4.950, Val: 27.040
Time: 4.960, Val: 27.296
Time: 4.970, Val: 27.554
Time: 4.980, Val: 27.813
Time: 4.990, Val: 28.072
Time: 5.000, Val: 28.333
|
Note that I have not checked the documentation. The approval is for the C++ code. |
|
Thank you Dr. Barra and @jbrzensk for checking my PR. The figure is now showing in the documentation. I also added the documentation to the index file as recommended. |
aboada
left a comment
aboada
left a comment
There was a problem hiding this comment.
Thank you for the changes, left some minor suggestions and an open question. Let me know if you have any questions!
examples/cpp/backward_euler.cpp
Outdated
| arma::mat y = arma::mat(1, t.size(), arma::fill::zeros); | ||
| constexpr double tol = 1e-6; // tolerance for fixed-point iteration | ||
|
|
||
There was a problem hiding this comment.
nit: I suggest removing introduced extra space. For formatting, I recommend clang-format though I'm not sure if that's well accepted within the library.
examples/cpp/backward_euler.cpp
Outdated
| break; | ||
| } | ||
| xn = xnp1; | ||
| // std::cout << k <<std::endl; // Print number of iterations until convergence |
There was a problem hiding this comment.
nit: I suggest removing unneeded commented line.
examples/cpp/backward_euler.cpp
Outdated
| if (std::abs(y(i+1) - y_old) < tol) { // Stopping criteria for approximate relative error | ||
| double xn = y_old; // Initial input for fixed-point iteration | ||
| double xnp1; // x at iteration n+1 | ||
| for (int n = 0; n < 100; n++) { // fixed-point iteration for rootfinding |
There was a problem hiding this comment.
I think the readability/maintainability of the code can improve by using a constant for the number of iterations (100). The constant can be declared as constexpr following the example of other constants in this file.
examples/cpp/backward_euler.cpp
Outdated
| double xnp1; // x at iteration n+1 | ||
| for (int n = 0; n < 100; 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 |
There was a problem hiding this comment.
What happens if we don't achieve the desired tolerance? Is that scenario handled in the current implementation? Should we inform the user?
There was a problem hiding this comment.
Hi Angel. The tolerance would be problem dependent (on the scale of the numbers in their problem) and needs to be specified by the user. If the tolerance isn't met, we can reduce the step-size of h or increase the tolerance. It's currently not handled in this implementation. However, I'll add a message to the user for situations where it does not converge.
examples/cpp/backward_euler.cpp
Outdated
| 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 iter = 100; // Max iterations for fixed-point iteration |
There was a problem hiding this comment.
nit: very minor but maybe max_iter can convey better the purpose of this variable.
aboada
left a comment
aboada
left a comment
There was a problem hiding this comment.
Thank you for applying clang-format! I noticed the comments at the top of the file were reformatted, likely because they were not marked as Doxygen comments (//!).
|
@valeriabarra Does the image integration work? I get an error ( missing image ). Do we need to add to the Makefile in the doc/sphinx in this section like we did for the MOLE OSE images? |
@jbrzensk , I think so, yes. Although now that I look more closely, it does not look like the approach we took in the Makefile is actually robust and scalable if we need to add more images to other sections of the docs (likely, in the Mathematical Framework section, for instance). If I look at libCEED's |
Off the top of my head, MyST's I've just published a PR (#260) that fixes this generically using Sphinx event hooks. It automatically collects images from multiple locations and fixes paths at both source and doctree levels. Originally I planned this for our long overdue of main branch RTD image invisibility, but it should serve dual purpose. I've tested it well on existing images, I'll keep in an eye on this PR to make sure it picks the untested images. |
5 participants
What type of PR is this? (check all applicable)
Description
I noticed that the backward Euler method was not implemented correctly. It is currently solving the ODE explicitly. The y at each iteration should be updated after each iteration instead. I also corrected some typos in the documentation.
Related Issues & Documents
examples/cpp/backward_euler.cpp
doc/sphinx/source/examples/Time-Integrators/backward_euler.md
QA Instructions, Screenshots, Recordings
Please replace this line with instructions on how to test your changes, a note
on the devices and browsers this has been tested on, as well as any relevant
images for UI changes.
Added/updated tests?
_We encourage you to test all code included with MOLE, including examples.
have not been included
Read Contributing Guide and Code of Conduct
[optional] Are there any post deployment tasks we need to perform?
[optional] What gif best describes this PR or how it makes you feel?