Skip to content

Commit

Permalink
WIP #1064 implement for determining default numerical integral method
Browse files Browse the repository at this point in the history
- If the expression has `Abs` function, use the `LegendreGauss` method

  Example input: `Integrate[Abs(x^2-2x), {x, -10, 10}] // N`
  Result should be `669.3282335875249`

- If the expression contains a variable that occurs in the exponent of
`Power` function, use the `GaussKronrod` method
  Example input: `x^x`, `3^(2x)`, `E^(-Sin(t))`

- Otherwise, use the `Romberg` method, since it is the optimized version
of trapezoid and Simpson methods
  • Loading branch information
axkr committed Sep 14, 2024
1 parent eac06e1 commit 2f35545
Show file tree
Hide file tree
Showing 2 changed files with 32 additions and 13 deletions.
Original file line number Diff line number Diff line change
Expand Up @@ -229,7 +229,7 @@ public NIntegrate() {
*/
@Override
public IExpr evaluate(final IAST ast, EvalEngine engine) {
String method = "LegendreGauss";
String method = "";
int maxPoints = DEFAULT_MAX_POINTS;
int maxIterations = DEFAULT_MAX_ITERATIONS;
int precisionGoal = 16; // automatic scale value
Expand Down Expand Up @@ -270,13 +270,23 @@ public IExpr evaluate(final IAST ast, EvalEngine engine) {
IAST list = (IAST) ast.arg2();
IExpr function = ast.arg1();
if (list.isAST3() && list.arg1().isSymbol()) {
IExpr x = list.arg1();
double min = list.arg2().evalf();
double max = list.arg3().evalf();
// if (min != null && max != null) {
if (function.isEqual()) {
IAST equalAST = (IAST) function;
function = F.Plus(equalAST.arg1(), F.Negate(equalAST.arg2()));
}
if (method.isEmpty()) {
method = "Romberg";
if (list.arg2().isDirectedInfinity() || list.arg3().isDirectedInfinity()) {
method = "LegendreGauss";
} else if (!function.isFree(a -> a == S.Abs || a == S.RealAbs, true)) {
method = "LegendreGauss";
} else if (!function.isFree(a -> a.isPower() && a.exponent().isFree(x), false)) {
method = "GaussKronrod";
}
}
try {
double result = integrate(method, list, min, max, function, maxPoints, maxIterations);
result = Precision.round(result, precisionGoal);
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -262,10 +262,10 @@ public void testIntegrate() {
// see github #120
check("Integrate(Ln(x)^2, {x,0,2})", //
"4-4*Log(2)+2*Log(2)^2");
check("Integrate(Ln(x)^2, {x,0,2}) // N", //
"2.18832");
check("NIntegrate(Ln(x)^2, {x,0,2}) // N", //
"2.1857");
checkNumeric("Integrate(Ln(x)^2, {x,0,2}) // N", //
"2.1883173055966214");
checkNumeric("NIntegrate(Ln(x)^2, {x,0,2}) // N", //
"2.188317305596199");

// see github #116
// should give (2*ArcTan((1 + 2*x)/Sqrt(3)))/Sqrt(3)
Expand Down Expand Up @@ -390,16 +390,21 @@ public void testNIntegrate() {
checkNumeric("NIntegrate(Cos(x), {x, 0, Pi})", //
"1.0E-16");
checkNumeric("NIntegrate(1/Sin(Sqrt(x)), {x, 0, 1}, PrecisionGoal->10)", //
"2.1108620052");
"2.1195255867");
}

@Test
public void testSystem171() {
check("N(Integrate(Sin(x),x))", "-Cos(x)");
check("N(Sin(x))", "Sin(x)");
check("Cancel((1*x+1/2*2)^((-1)*2)*1^(-1)^(-1))", "1/(1+x)^2");
check("Integrate(1/(a+b*x),x)", "Log(a+b*x)/b");
check("Integrate((a+b*x)^(1/3),x)", "3/4*(a+b*x)^(4/3)/b");
check("N(Integrate(Sin(x),x))", //
"-Cos(x)");
check("N(Sin(x))", //
"Sin(x)");
check("Cancel((1*x+1/2*2)^((-1)*2)*1^(-1)^(-1))", //
"1/(1+x)^2");
check("Integrate(1/(a+b*x),x)", //
"Log(a+b*x)/b");
check("Integrate((a+b*x)^(1/3),x)", //
"3/4*(a+b*x)^(4/3)/b");
}

@Test
Expand Down Expand Up @@ -610,11 +615,15 @@ public void testXReciprocalIssue1064() {
check("NIntegrate(1/x, {x,0,1},Method->\"LegendreGauss\")", //
"10.37476");
check("N(Integrate(1/x, {x,0,1}))", //
"10.37476");
"Integrate(1/x,{x,0,1})");

// message - Integrate: Integral of 1/x does not converge on {x,0,1}.
check("Integrate(1/x, {x,0,1})", //
"Integrate(1/x,{x,0,1})");

checkNumeric("Integrate(Abs(x^2-2*x), {x, -10, 10}) // N", //
"669.3282335875249");

}

}

0 comments on commit 2f35545

Please sign in to comment.