]>
Commit | Line | Data |
---|---|---|
92f5a8d4 TL |
1 | // Copyright Matthew Pulver 2018 - 2019. |
2 | // Distributed under the Boost Software License, Version 1.0. | |
3 | // (See accompanying file LICENSE_1_0.txt or copy at | |
4 | // https://www.boost.org/LICENSE_1_0.txt) | |
5 | ||
6 | #include <boost/math/differentiation/autodiff.hpp> | |
7 | #include <boost/multiprecision/cpp_bin_float.hpp> | |
8 | #include <iostream> | |
9 | ||
10 | using namespace boost::math::differentiation; | |
11 | ||
12 | template <typename W, typename X, typename Y, typename Z> | |
13 | promote<W, X, Y, Z> f(const W& w, const X& x, const Y& y, const Z& z) { | |
14 | using namespace std; | |
15 | return exp(w * sin(x * log(y) / z) + sqrt(w * z / (x * y))) + w * w / tan(z); | |
16 | } | |
17 | ||
18 | int main() { | |
19 | using float50 = boost::multiprecision::cpp_bin_float_50; | |
20 | ||
21 | constexpr unsigned Nw = 3; // Max order of derivative to calculate for w | |
22 | constexpr unsigned Nx = 2; // Max order of derivative to calculate for x | |
23 | constexpr unsigned Ny = 4; // Max order of derivative to calculate for y | |
24 | constexpr unsigned Nz = 3; // Max order of derivative to calculate for z | |
25 | // Declare 4 independent variables together into a std::tuple. | |
26 | auto const variables = make_ftuple<float50, Nw, Nx, Ny, Nz>(11, 12, 13, 14); | |
27 | auto const& w = std::get<0>(variables); // Up to Nw derivatives at w=11 | |
28 | auto const& x = std::get<1>(variables); // Up to Nx derivatives at x=12 | |
29 | auto const& y = std::get<2>(variables); // Up to Ny derivatives at y=13 | |
30 | auto const& z = std::get<3>(variables); // Up to Nz derivatives at z=14 | |
31 | auto const v = f(w, x, y, z); | |
32 | // Calculated from Mathematica symbolic differentiation. | |
33 | float50 const answer("1976.319600747797717779881875290418720908121189218755"); | |
34 | std::cout << std::setprecision(std::numeric_limits<float50>::digits10) | |
35 | << "mathematica : " << answer << '\n' | |
36 | << "autodiff : " << v.derivative(Nw, Nx, Ny, Nz) << '\n' | |
37 | << std::setprecision(3) | |
38 | << "relative error: " << (v.derivative(Nw, Nx, Ny, Nz) / answer - 1) << '\n'; | |
39 | return 0; | |
40 | } | |
41 | /* | |
42 | Output: | |
43 | mathematica : 1976.3196007477977177798818752904187209081211892188 | |
44 | autodiff : 1976.3196007477977177798818752904187209081211892188 | |
45 | relative error: 2.67e-50 | |
46 | **/ |