diff options
| author | Tavian Barnes <tavianator@gmail.com> | 2010-10-11 02:51:41 -0400 |
|---|---|---|
| committer | Tavian Barnes <tavianator@gmail.com> | 2010-10-11 02:51:41 -0400 |
| commit | f123be63a5760782c563d81d7842de6cc3e5646b (patch) | |
| tree | ceb533b40e917156cad556ca99c75ebe65ae2d21 | |
| parent | 2c2b69e1df183f118c45d74c18c7e84934aaff1f (diff) | |
| download | vz-f123be63a5760782c563d81d7842de6cc3e5646b.tar.xz | |
Support adaptive integration with equation systems.
| -rw-r--r-- | src/vZ/EquationSystem.hpp | 14 | ||||
| -rw-r--r-- | tests/EquationSystem.cpp | 10 |
2 files changed, 20 insertions, 4 deletions
diff --git a/src/vZ/EquationSystem.hpp b/src/vZ/EquationSystem.hpp index 0963baf..6409f18 100644 --- a/src/vZ/EquationSystem.hpp +++ b/src/vZ/EquationSystem.hpp @@ -21,6 +21,8 @@ #ifndef VZ_EQUATIONSYSTEM_HPP #define VZ_EQUATIONSYSTEM_HPP +#include <algorithm> +#include <cmath> #include <cstddef> namespace vZ @@ -149,6 +151,18 @@ namespace vZ } return *this; } + + template <std::size_t N, typename T> + typename EquationSystem<N, T>::Scalar + abs(const EquationSystem<N, T>& es) + { + typename EquationSystem<N, T>::Scalar ret(0); + for (std::size_t i = 0; i < N; ++i) { + using std::abs; + ret = std::max(ret, abs(es[i])); + } + return ret; + } } #endif // VZ_EQUATIONSYSTEM_HPP diff --git a/tests/EquationSystem.cpp b/tests/EquationSystem.cpp index 973dd30..497a065 100644 --- a/tests/EquationSystem.cpp +++ b/tests/EquationSystem.cpp @@ -22,8 +22,8 @@ main() Y y; y[0] = 1.0; y[1] = 1.0; - vZ::GenericEulerIntegrator<Y> integrator(f); - integrator.y(y).x(0.0).h(0.01); + vZ::GenericDP45Integrator<Y> integrator(f); + integrator.tol(1e-6).y(y).x(0.0).h(0.06); integrator.integrate(2.0); @@ -33,10 +33,12 @@ main() std::cout << std::setprecision(10) << "Numerical: " << actual << std::endl << "Expected: " << expected << std::endl - << "Iterations: " << integrator.iterations() << std::endl; + << "h: " << integrator.h() << std::endl + << "Iterations: " << integrator.iterations() << std::endl + << "Rejections: " << integrator.rejections() << std::endl; double error = std::fabs(expected - actual)/expected; - if (error > 0.01) { + if (error > 6.0e-7) { std::cerr << "Error: " << 100.0*error << "%" << std::endl; return EXIT_FAILURE; } else { |