Взаимодействия среды

Функция чувствительна к переменной окружения DIGITS, который определяет числовую рабочую точность.

Пример 1

Мы демонстрируем стандартные призывы к адаптивному численному интегрированию:

numeric::quadrature(exp(x^2), x = -1..1)

numeric::quadrature(max(1/10, cos(PI*x)), x = -2..0.0123)

numeric::quadrature(exp(-x^2), x = -2..infinity)

Цель точности устанавливается DIGITS:

DIGITS := 50:
numeric::quadrature(besselJ(0, x), x = 0..PI)

Обратите внимание на то, что из-за внутреннего адаптивного механизма, выбор различных методов не должен иметь никакого значительного эффекта на квадратурный результат:

DIGITS := 10:
numeric::quadrature(sin(x)/x, x = -1..10, GaussLegendre = 5),
numeric::quadrature(sin(x)/x, x = -1..10, GaussLegendre = 160),
numeric::quadrature(sin(x)/x, x = -1..10, NewtonCotes = 8)

Пример 2

Пользователь должен убедиться, что подынтегральное выражение четко определено и достаточно регулярное. Следующие сбои, потому что квадратура Коутса Ньютона пытается оценить подынтегральное выражение на контурах:

numeric::quadrature(sin(x)/x, x = 0..1, NewtonCotes = 8)

Error: Division by zero. [_power]
  Evaluating: Quadsum

Можно исправить эту проблему присвоить значение f(0). Подынтегральное выражение передается интегратору как hold(f), чтобы предотвратить преждевременную оценку f(x) к sin(x)/x. Внутренне, numeric::quadrature заменяет x численными значениями и затем оценивает подынтегральное выражение. Обратите внимание на то, что нужно задать f(0.0) := 1, а не f(0) := 1:

f := x -> sin(x)/x:
f(0.0) := 1:
numeric::quadrature(hold(f)(x), x = 0..1, NewtonCotes = 8)

Метод по умолчанию (Квадратура Гаусса - Лежандра) не оценивает подынтегральное выражение в конечных точках:

numeric::quadrature(sin(x)/x, x = 0..1)

Тем не менее, проблемы могут все еще возникнуть для несобственных интегралов:

numeric::quadrature(ln((1 + x^4)^2 - 1), x = 0 .. 1)

Warning: Precision goal not achieved after 10000 function calls. Increase 'MaxCalls' and try again for a more accurate result. [numeric::quadrature]

В этом примере подынтегральное выражение оценено близко к 0. Выражение (1 + x ⁴) ² - 1 страдает от серьезной числовой отмены и во власти округления. Численно стабильное представление (1 + x ⁴) ² - 1 = x ⁴   (x ⁴ + 2):

numeric::quadrature(ln(x^4*(x^4 + 2)), x = 0..1)

Обратите внимание на то, что сходимость является довольно медленной, потому что подынтегральное выражение не ограничено.

delete f:

Пример 3

Мы демонстрируем многомерную квадратуру:

Q := numeric::quadrature:
Q(Q(exp(x*y), x = 0..y), y = 0..1)

Также более составные типы вложенных вызовов возможны. Мы используем численно заданные функции

b := y -> Q(exp(-t^2-t*y), t = y..infinity):

и

f := y -> cos(y^2) +  Q(exp(x*y), x = 0..b(y)):

для следующей квадратуры:

Q(f(y), y = 0..1)

Много размерная квадратура является в вычислительном отношении дорогой. Обратите внимание на то, что вызов numeric::quadrature оценивает подынтегральное выражение, по крайней мере, 3 n времена, где n является количеством узлов внутреннего квадратурного правила (по умолчанию, n = 20 для DIGITS ≤ 10). Следующая тройная квадратура вызвала бы функцию exp (не менее чем 3 20) ³ = 216000 раз!

Q(Q(Q(exp(x*y*z), x = 0..y+z), y = 0..z), z = 0..1)

Для низких целей точности достаточны квадратурные правила низкоуровневые. В следующем мы уменьшаем вычислительные затраты при помощи квадратуры Гаусса - Лежандра с 5 узлами. Мы используем краткое обозначение GL, чтобы задать метод GaussLegendre:

DIGITS := 4: 
Q(Q(Q(exp(x*y*z), x=0..y+z, GL=5), y=0..z, GL=5), z=0..1, GL=5)

delete Q, b, f, DIGITS:

Пример 4

Мы демонстрируем, как должны быть обработаны подынтегральные выражения, данные пользовательскими процедурами. Следующее подынтегральное выражение

f := proc(x) begin 
      if x<1 then sin(x^2) else cos(x^5) end_if
    end_proc:

не может быть вызван символьным аргументом:

f(x)

Error: Unable to evaluate to Boolean. [_less]
  Evaluating: f

Следовательно, нужно использовать hold, чтобы предотвратить преждевременную оценку f(x):

numeric::quadrature(hold(f)(x), x = -1..PI/2)

Обратите внимание на то, что вышеупомянутое подынтегральное выражение прерывисто в x = 1, вызывая медленную сходимость числовой квадратуры. Намного более эффективно разделить интеграл в две подквадратуры со сглаженными подынтегральными выражениями:

numeric::quadrature(sin(x^2), x = -1..1) +
numeric::quadrature(cos(x^5), x = 1..PI/2)

Смотрите Пример 5 для многомерной квадратуры пользовательских процедур.

delete f:

Пример 5

Следующее подынтегральное выражение

f := proc(x, y) begin 
       if x<y then x-y else (x-y) + (x-y)^5 end_if
     end_proc:

может только быть вызван числовыми аргументами и должен быть задержан дважды для 2-мерной квадратуры:

Q := numeric::quadrature:

Q(Q(hold(hold(f))(x, y), x = 0..1), y = 0..1)

Обратите внимание на то, что задержка подынтегрального выражения через hold не будет работать на тройную или более многомерную квадратуру! Пользователь может обработать это путем проверки, что подынтегральное выражение может также быть оценено для символьных аргументов:

f := proc(x, y, z)
     begin
       if not testtype([args()], Type::ListOf(Type::Numeric))
          then return(procname(args()))
       end_if;
       if x^2 + y^2 + z^2 <= 1 
          then return(1) 
          else return(0) 
       end_if
     end_proc:

Обратите внимание на то, что эта функция не непрерывна, подразумевая медленную сходимость числовой квадратуры. Поэтому мы используем низкую цель точности только 2 цифр и уменьшаем затраты при помощи квадратурного правила низкоуровневого:

DIGITS := 2: 
Q(Q(Q(f(x, y, z), x=0..1, GL=5), y=0..1, GL=5), z=0..1, GL=5)

delete f, Q, DIGITS:

Пример 6

Следующий пример использует неадаптивную квадратуру Гаусса-Tschebyscheff с растущим числом узлов. Результаты сходятся быстро к точному значению:

a := exp(x)/sqrt(1 - x^2), x = -1..1:
numeric::quadrature(a, Adaptive = FALSE, GT = n) $ n = 3..7

delete a:

Пример 7

Несобственный интеграл существует. Числовая сходимость, однако, является довольно медленной из-за особенности в x = 0:

numeric::quadrature(x^(-9/10), x = 0..1)

Warning: Precision goal not achieved after 10000 function calls. Increase 'MaxCalls' and try again for a more accurate result. [numeric::quadrature]

Мы удаляем предел для количества вызовов функции и позволяем numeric::quadrature сточиться вперед, пока результат не найден. Время для вычисления растет соответственно, последняя цифра является неправильной должная округлить эффекты:

numeric::quadrature(x^(-9/10), x = 0..1, MaxCalls = infinity)

Пример 8

Следующий интеграл не существует в Римановом смысле, потому что подынтегральное выражение не ограничено:

numeric::quadrature(1/x, x = -1..2)

Warning: Precision goal not achieved after 10000 function calls. Increase 'MaxCalls' and try again for a more accurate result. [numeric::quadrature]

Мы увеличиваем MaxCalls. Нет никакой сходимости числового алгоритма, потому что интеграл не существует. Следовательно, некоторая внутренняя проблема должна возникнуть: numeric::quadrature достигает своей максимальной рекурсивной глубины:

numeric::quadrature(1/x, x = -1..2, MaxCalls = infinity)

Warning: Precision goal not achieved after 'MAXDEPTH=500' recursive calls. There might be a singularity of '1/x' close to 'x=3.910318545e-148'. Increase 'MAXDEPTH' and try again for a more accurate result. [adaptiveQuad]

Возвращаемые значения

Число с плавающей точкой возвращено, если нечисловые символьные объекты в подынтегральном выражении f(x) или в контурах a,b не предотвращают численную оценку. В этом случае numeric::quadrature возвращает себя символически. Если числовые проблемы происходят, то FAIL возвращен.