diff --git a/include/expolit/differentiation.hh b/include/expolit/differentiation.hh index 5ad6b81..2bb872f 100644 --- a/include/expolit/differentiation.hh +++ b/include/expolit/differentiation.hh @@ -1,87 +1,123 @@ /** * @file * * @author Lucas Frérot * * @section LICENSE * * Copyright (©) 2018 EPFL (Ecole Polytechnique Fédérale de * Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des * Solides) * * Expolit is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Expolit is distributed in the hope that it will be useful, but WITHOUT * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS * FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for * more details. * * You should have received a copy of the GNU Lesser General Public License * along with Expolit. If not, see . * */ #ifndef DIFFERENTIATION_HH #define DIFFERENTIATION_HH #include "exponential.hh" +#include "litteral.hh" #include "operations.hh" #include "polynomial.hh" -#include "litteral.hh" #include "types.hh" #include namespace expolit { /// Differentiate polynomial template constexpr auto differentiate(const Polynomial& p) { Polynomial dp; for (UInt i = 0; i < dp.terms; ++i) dp.coeffs[i] = p.coeffs[i + 1] * static_cast(i + 1); return dp; } /// Differentiate constant template constexpr auto differentiate(const Constant& p) { return Constant({0}); } /// Differentiate litteral template constexpr auto differentiate(const Litteral& p) { return 0; } /// Differentiate exponential template constexpr auto differentiate(const Exponential& e) { return differentiate(e.expression) * e; } /// Differentiate generic product template constexpr auto differentiate(const Product& p) { auto&& o = p.operands; return differentiate(o.first) * o.second + o.first * differentiate(o.second); } /// Differentiate constant + expression template constexpr auto differentiate(const Sum, U>& p) { return differentiate(p.operands.second); } +/// Differentiate expression + constant +template +constexpr auto differentiate(const Sum>& p) { + return differentiate(p.commute()); +} + +/// Differentiate litteral + expression +template +constexpr auto differentiate(const Sum, U>& p) { + return differentiate(p.operands.second); +} + +/// Differentiate expression + litteral +template +constexpr auto differentiate(const Sum>& p) { + return differentiate(p.commute()); +} + /// Differentiate constant * expression template constexpr auto differentiate(const Product, U>& p) { return p.operands.first * differentiate(p.operands.second); } +/// Differentiate expression * constant +template +constexpr auto differentiate(const Product>& p) { + return differentiate(p.commute()); +} + +/// Differentiate litteral * expression +template +constexpr auto differentiate(const Product, U>& p) { + return p.operands.first * differentiate(p.operands.second); +} + +/// Differentiate expression * litteral +template +constexpr auto differentiate(const Product>& p) { + return differentiate(p.commute()); +} + } // namespace expolit #endif diff --git a/tests/tests_gtest.cc b/tests/tests_gtest.cc index 2453060..c81c10b 100644 --- a/tests/tests_gtest.cc +++ b/tests/tests_gtest.cc @@ -1,156 +1,162 @@ /** * @file * * @author Lucas Frérot * * @section LICENSE * * Copyright (©) 2018 EPFL (Ecole Polytechnique Fédérale de * Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des * Solides) * * Expolit is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Expolit is distributed in the hope that it will be useful, but WITHOUT * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS * FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for * more details. * * You should have received a copy of the GNU Lesser General Public License * along with Expolit. If not, see . * */ #include #include #include using namespace expolit; static constexpr Polynomial Zi({0, 1}); static constexpr Polynomial, 1> Zc({0, 1}); static constexpr Polynomial Zr({0, 1}); using lin = std::remove_cv_t; template constexpr void compare(const Polynomial& p1, const Polynomial& p2) { for (UInt i = 0; i < p1.terms; ++i) ASSERT_NEAR(std::abs(p1.coeffs[i] - p2.coeffs[i]), 0, 1e-14) << "Polynomials are not equal (i = " << i << ")"; } template constexpr void compare(const Polynomial& p1, const Polynomial& p2) { for (UInt i = 0; i < p1.terms; ++i) ASSERT_EQ(p1.coeffs[i], p2.coeffs[i]) << "Polynomials are not equal (i = " << i << ")"; } TEST(Polynomial, initialization) { constexpr auto p = 1 + 2 * pow<2>(Zi) + pow<3>(Zi); constexpr Int z = 3; constexpr Int sol = 1 + 2 * z * z + z * z * z; constexpr Int res = p(z); static_assert(sol == res, "Bad coefficients"); // This guy cannot be consexpr because std::complex is stupid auto p2 = 1. + Zc * Zc; std::complex z2(0, 1); ASSERT_NEAR(std::abs(p2(z2)), 0, 1e-14) << "Bad coefficents (complex)"; } TEST(Polynomial, sum) { constexpr auto p = 1 + 2 * pow<2>(Zi) + pow<3>(Zi); constexpr auto p2 = Zi + 2 * Zi * Zi; constexpr auto res = 1 + Zi + 4 * Zi * Zi + Zi * Zi * Zi; compare(p + p2, res); } TEST(Polynomial, product) { constexpr auto p = 1 + 2 * pow<2>(Zi) + pow<3>(Zi); constexpr auto p2 = Zi + 2 * pow<2>(Zi); constexpr auto res = Zi + 2 * pow<2>(Zi) + 2 * pow<3>(Zi) + 5 * pow<4>(Zi) + 2 * pow<5>(Zi); compare(p * p2, res); } TEST(Polynomial, derivative) { constexpr auto p = 1 + 2 * pow<2>(Zi) + pow<3>(Zi); constexpr auto res = 4 * Zi + 3 * Zi * Zi; compare(differentiate(p), res); } TEST(Polynomial, integration) { constexpr auto p = 1. + 2. * pow<2>(Zr) + pow<3>(Zr); constexpr auto res = Zr + 2 / 3. * pow<3>(Zr) + 1 / 4. * pow<4>(Zr); compare(integrate(p), res); } TEST(Exponential, evaluation) { constexpr auto e = exp(2. * Zr); constexpr Real z = 3.4; const Real res = e(z); // cannot be constexpr because of exponential ASSERT_NEAR(std::exp(2 * z), res, 1e-14) << "Exponential evaluation fail"; } TEST(Exponential, derivative) { constexpr auto e = exp(2. * Zr); constexpr auto diff = differentiate(e); constexpr auto res = Constant({2}) * e; compare(res.operands.first, diff.operands.first); compare(e.expression, diff.operands.second.expression); } TEST(Exponential, integration) { constexpr auto e = exp(2. * Zr); constexpr auto integ = integrate(e); constexpr auto res = Constant({1 / 2.}) * e; compare(res.operands.first, integ.operands.first); compare(e.expression, integ.operands.second.expression); } TEST(Integration, parts) { constexpr auto e = exp(2. * Zr); constexpr auto expr = e * (2. * Zr); constexpr auto res = Zr * e - Constant({1 / 2.}) * e; constexpr auto integ = integrate(expr); compare(integ.operands.first.operands.first, res.operands.first.operands.first); compare(integ.operands.second.operands.first, res.operands.second.operands.first); } TEST(Integration, definite_integral) { constexpr auto p = Zr; constexpr std::pair bounds(0, 1); constexpr auto res = definite_integral(bounds, p); ASSERT_NEAR(0.5, res, 1e-14) << "Definite integral fail"; } using Q = Litteral; TEST(Litteral, substitution) { constexpr Q q; constexpr auto p = q * Zi; constexpr auto c = Constant({2}); constexpr auto p_sub = substitute(c, p); compare(2 * Zi, p_sub); constexpr auto expo = exp(q * Zr); constexpr auto x_sub = substitute(c, expo); compare(2. * Zr, x_sub.expression); } TEST(Litteral, differentiation) { - + constexpr Q q; + constexpr auto c = Constant({2}); + constexpr auto expo = exp(q * Zr); + constexpr auto dexpo = differentiate(expo); + constexpr auto res = substitute(c, dexpo); + compare(c, res.operands.first); + compare(2. * Zr, res.operands.second.expression); } TEST(Litteral, integration) { }