diff --git a/src/Exponent.cpp b/src/Exponent.cpp index 08e3bd54..65e624d2 100644 --- a/src/Exponent.cpp +++ b/src/Exponent.cpp @@ -51,6 +51,18 @@ auto Exponent::Integrate(const Expression& integrationVariable) cons return adder.Accept(simplifyVisitor).value(); } } + + // EulerNumber raised to a (degree-one) Variable + if (auto eulerBase = RecursiveCast>(*simplifiedExponent); eulerBase != nullptr) { + const Variable& expPow = eulerBase->GetLeastSigOp(); + + Add adder { + Exponent { EulerNumber {}, expPow }, + Variable { "C" } + }; + + return adder.Accept(simplifyVisitor).value(); + } } Integral integral { *(this->Copy()), *(integrationVariable.Copy()) }; diff --git a/src/Multiply.cpp b/src/Multiply.cpp index cbebebc3..dd522474 100644 --- a/src/Multiply.cpp +++ b/src/Multiply.cpp @@ -5,11 +5,14 @@ #include "Oasis/Multiply.hpp" #include "Oasis/Add.hpp" #include "Oasis/Divide.hpp" +#include "Oasis/EulerNumber.hpp" #include "Oasis/Exponent.hpp" #include "Oasis/Imaginary.hpp" #include "Oasis/Integral.hpp" +#include "Oasis/Log.hpp" #include "Oasis/Matrix.hpp" #include "Oasis/RecursiveCast.hpp" +#include "Oasis/Subtract.hpp" #define EPSILON 10E-6 @@ -42,7 +45,101 @@ auto Multiply::Integrate(const Expression& integrationVariable) cons } } - // TODO: Implement integration by parts + // Apply Integration By Parts + else { + // Start by assuming u is MostSigOp and dv is LeastSigOp + auto multiplyCopy = std::make_unique>(Multiply { this->GetMostSigOp(), this->GetLeastSigOp() }); + + // Check by LIPET if it's necessary to swap operands so u is LeastSigOp and dv is MostSigOp + // LIPET: Logarithm case + if (!(multiplyCopy->GetMostSigOp().Is>()) + && multiplyCopy->GetLeastSigOp().Is>()) { + multiplyCopy = std::make_unique>(multiplyCopy->SwapOperands()); + } + + // LIPET: Inverse trigonometry is not implemented yet in Oasis + + // LIPET: Swap polynomial (linear variable) and EulerNumber + if ((multiplyCopy->GetLeastSigOp().Is() + && multiplyCopy->GetMostSigOp().Is>())) { + multiplyCopy = std::make_unique>(multiplyCopy->SwapOperands()); + } + + // LIPET: Swap polynomial (exponential) and EulerNumber + if (auto polynomial = RecursiveCast>(multiplyCopy->GetLeastSigOp()); polynomial != nullptr) { + if (auto euler = RecursiveCast>(multiplyCopy->GetMostSigOp()); euler != nullptr) { + multiplyCopy = std::make_unique>(multiplyCopy->SwapOperands()); + } + } + + // LIPET: Trigonometry is not implemented yet in Oasis + + // Integrate dv and differentiate u to attain v and du + auto v = multiplyCopy->GetLeastSigOp().Integrate(integrationVariable); + auto du = multiplyCopy->GetMostSigOp().Differentiate(integrationVariable); + + // Remove the +C from v + auto vPlusC = RecursiveCast>(*v); + v = vPlusC->GetMostSigOp().Copy(); + + // Initialize the coefficients of v and du to 1 + auto vCoefficient = std::make_unique(Real { 1 }); + auto duCoefficient = std::make_unique(Real { 1 }); + + // Make a copy of v so that the coefficient can be factored out when computing vdu, + // while preserving the original v for the multiplication of u*v + auto v_VDU = v->Copy(); + + // Factor out the coefficient for v; v_VDU now has a coefficient of one + if (auto coefficientMultiply = RecursiveCast>(*v); coefficientMultiply != nullptr) { + vCoefficient = RecursiveCast(coefficientMultiply->GetMostSigOp()); + v_VDU = coefficientMultiply->GetLeastSigOp().Accept(simplifyVisitor).value(); + } + + // Factor out the coefficient for du; du now has a coefficient of one + if (auto coefficientMultiply = RecursiveCast>(*du); coefficientMultiply != nullptr) { + duCoefficient = RecursiveCast(coefficientMultiply->GetMostSigOp()); + du = coefficientMultiply->GetLeastSigOp().Accept(simplifyVisitor).value(); + } + + // Attain a constant for the coefficient of vdu + auto vduCoefficient = std::make_unique>(Multiply { *vCoefficient, *duCoefficient })->Accept(simplifyVisitor).value(); + + // Multiply v and du to attain partVDU, which is vdu with a coefficient of one + auto partVDU = std::make_unique>(Multiply { *v_VDU, *du })->Accept(simplifyVisitor).value(); + + // Prevent infinite recursion + // If fullVDU is equal to the original Multiply instance, + // then IBP cannot reduce the integrand, so return an integral of the Multiply instance + auto fullVDU = std::make_unique>(Multiply { *vduCoefficient, *partVDU }); + if (multiplyCopy->Equals(*fullVDU)) { + Integral integral { *(multiplyCopy->Copy()), *(integrationVariable.Copy()) }; + return integral.Copy(); + } + + // If partVDU is a Multiply instance it will use IBP again, + // otherwise it will find another integration method to use. + // Implicit base case: IBP will eventually reduce the integrand such that it is not a Multiply instance. + auto integratedPartVDU = partVDU->Integrate(integrationVariable); + + // Remove the +C from the integration result + if (auto removeC = RecursiveCast>(*integratedPartVDU); removeC != nullptr) { + integratedPartVDU = removeC->GetMostSigOp().Copy(); + } + + // Apply coefficient correction for integratedPartVDU, thus building integratedFullVDU + auto integratedFullVDU = std::make_unique>(Multiply { *vduCoefficient, *integratedPartVDU }); + + // Apply the Integration By Parts formula, and add the +C at the end + Add, Expression>, Variable> adder { + Subtract, Expression> { + Multiply { multiplyCopy->GetMostSigOp(), *v }, + *integratedFullVDU }, + Variable { "C" } + }; + + return adder.Accept(simplifyVisitor).value(); + } } Integral integral { *(this->Copy()), *(integrationVariable.Copy()) }; diff --git a/tests/IntegrateTests.cpp b/tests/IntegrateTests.cpp index 7e1c2ede..4ef48f57 100644 --- a/tests/IntegrateTests.cpp +++ b/tests/IntegrateTests.cpp @@ -2,7 +2,9 @@ #include "Oasis/Add.hpp" #include "Oasis/Divide.hpp" +#include "Oasis/EulerNumber.hpp" #include "Oasis/Exponent.hpp" +#include "Oasis/Log.hpp" #include "Oasis/Multiply.hpp" #include "Oasis/Real.hpp" #include "Oasis/Subtract.hpp" @@ -181,4 +183,160 @@ TEST_CASE("Integrate Add Rule Like Terms", "[Integrate][Add][Like]") auto simplified = integrated->Accept(simplifyVisitor).value(); REQUIRE(simplified->Equals(*(integral.Accept(simplifyVisitor).value()))); +} + +TEST_CASE("Integrate EulerNumber Raised to Power x", "[Integrate][Euler]") +{ + Oasis::Variable var { "x" }; + + Oasis::Exponent integrand { Oasis::EulerNumber{}, Oasis::Variable { var.GetName() } }; + + Oasis::Add, Oasis::Variable> integral { + Oasis::Add { + Oasis::Exponent { Oasis::EulerNumber{}, Oasis::Variable { var.GetName() } }, + Oasis::Variable { "C" } + } + }; + + auto ptr = integral.Accept(simplifyVisitor).value(); + auto integrated = integrand.Integrate(var); + REQUIRE((integrated->Equals(*ptr))); +} + +TEST_CASE("Integration By Parts: Variable and Euler's Number", "[Integrate][Variable][Euler]") +{ + Oasis::Variable var { "x" }; + + Oasis::Multiply> integrand { + Oasis::Variable { var.GetName() }, + Oasis::Exponent { Oasis::EulerNumber{}, Oasis::Variable { var.GetName() } } + }; + + Oasis::Add>, Oasis::Exponent>, Oasis::Variable> integral { + Oasis::Add { + Oasis::Subtract { + Oasis::Multiply { + Oasis::Variable { var.GetName() }, + Oasis::Exponent { Oasis::EulerNumber{}, Oasis::Variable { var.GetName() } } }, + Oasis::Exponent { Oasis::EulerNumber{}, Oasis::Variable { var.GetName() } } }, + Oasis::Variable { "C" } + } + }; + + auto ptr = integral.Accept(simplifyVisitor).value(); + auto integrated = integrand.Integrate(var); + REQUIRE((integrated->Equals(*ptr))); + + integrated = integrand.SwapOperands().Integrate(var); + REQUIRE((integrated->Equals(*ptr))); +} + +TEST_CASE("Integration By Parts: Exponent and Euler's Number", "[Integrate][Exponent][Euler]") +{ + Oasis::Variable var { "x" }; + + Oasis::Multiply, Oasis::Exponent> integrand { + Oasis::Exponent { Oasis::Variable { var.GetName() }, Oasis::Real { 2 } }, + Oasis::Exponent { Oasis::EulerNumber{}, Oasis::Variable { var.GetName() } } + }; + + Oasis::Add, + Oasis::Add, Oasis::Multiply>>>, Oasis::Variable> integral { + Oasis::Add { + Oasis::Multiply { + Oasis::Exponent { Oasis::EulerNumber{}, Oasis::Variable { var.GetName() } }, + Oasis::Add { + Oasis::Exponent { Oasis::Variable { var.GetName() }, Oasis::Real { 2 } }, + Oasis::Multiply { + Oasis::Real {-2}, + Oasis::Add { + Oasis::Variable { var.GetName() }, + Oasis::Real { -1 } + } + }, + } + }, + Oasis::Variable { "C" } + } + }; + + + auto ptr = integral.Accept(simplifyVisitor).value(); + auto integrated = integrand.Integrate(var); + REQUIRE((integrated->Equals(*ptr))); + + integrated = integrand.SwapOperands().Integrate(var); + REQUIRE((integrated->Equals(*ptr))); + +} + +TEST_CASE("Integration By Parts: Variable and Logarithm", "[Integrate][Variable][Logarithm]") +{ + Oasis::Variable var { "x" }; + + Oasis::Multiply> integrand { + Oasis::Variable { var.GetName() }, + Oasis::Log { Oasis::EulerNumber{}, Oasis::Variable { var.GetName() } } + }; + + Oasis::Add, + Oasis::Add, + Oasis::Log>, Oasis::Divide>>, Oasis::Variable> integral { + Oasis::Add { + Oasis::Multiply { + Oasis::Exponent { Oasis::Variable { var.GetName() }, Oasis::Real { 2 } }, + Oasis::Add { + Oasis::Multiply { + Oasis::Divide { Oasis::Real { 1 }, Oasis::Real { 2 } }, + Oasis::Log { Oasis::EulerNumber{}, Oasis::Variable { var.GetName() } } + }, + Oasis::Divide { Oasis::Real { -1 }, Oasis::Real { 4 } }, + } + }, + Oasis::Variable { "C" } + } + }; + + auto ptr = integral.Accept(simplifyVisitor).value(); + auto integrated = integrand.Integrate(var); + REQUIRE((integrated->Equals(*ptr))); + + integrated = integrand.SwapOperands().Integrate(var); + REQUIRE((integrated->Equals(*ptr))); +} + +TEST_CASE("Integration By Parts: Exponent and Logarithm", "[Integrate][Exponent][Logarithm]") +{ + Oasis::Variable var { "x" }; + + Oasis::Multiply, Oasis::Log> integrand { + Oasis::Exponent { Oasis::Variable { var.GetName() }, Oasis::Real { 2 } }, + Oasis::Log { Oasis::EulerNumber{}, Oasis::Variable { var.GetName() } } + }; + + Oasis::Add, + Oasis::Add, + Oasis::Log>, Oasis::Divide>>, Oasis::Variable> integral { + Oasis::Add { + Oasis::Multiply { + Oasis::Exponent { Oasis::Variable { var.GetName() }, Oasis::Real { 3 } }, + Oasis::Add { + Oasis::Multiply { + Oasis::Divide { Oasis::Real { 1 }, Oasis::Real { 3 } }, + Oasis::Log { Oasis::EulerNumber{}, Oasis::Variable { var.GetName() } } + }, + Oasis::Divide { Oasis::Real { -1 }, Oasis::Real { 9 } }, + } + }, + Oasis::Variable { "C" } + } + }; + + auto ptr = integral.Accept(simplifyVisitor).value(); + auto integrated = integrand.Integrate(var); + REQUIRE((integrated->Equals(*ptr))); + + integrated = integrand.SwapOperands().Integrate(var); + REQUIRE((integrated->Equals(*ptr))); } \ No newline at end of file