Иллюстрированный самоучитель по Mathematica
Численные расчеты — пакет NumericalMath
Пакет расширения
NumericalMath содержит множество полезных функций для тех, кто имеет дело с
численными расчетами. В их числе функции для выполнения высокоточных аппроксимаций
рациональными функциями, численного интегрирования и дифференцирования, вычисления
пределов функций, решения уравнений, разложения в ряд и т. д. Ниже описано подавляющее
большинство функций этого расширения. Исключены лишь отдельные функции, представляющие
ограниченный интерес и несложные для самостоятельного изучения (в подпаке-mах
Butcher, Microscope и ComputerArithmetic).
Аппроксимация аналитических функций — Approximations
Подпакет
Approximations содержит ряд функций для улучшенной рациональной аппроксимации
аналитических функций. Для рациональной интерполяции и аппроксимации функций
по заданным значениям абсцисс служит следующая функция:
-
Rationallnterpolation
[f, {x,m, k}, {x
1
, x
2
, ...,.x
m+k+1
} ] — возвращает
аппроксимирующее функцию f выражение в виде отношения полиномов а степенью
полинома числителя m и знаменателя k в абсциссах, заданных списком {x
l
,x
2
,...,x
m+jt+1
}.
Пример применения
этой функции:
<<NumericalMath `Approximations`
ril =
Rationallnterpolation[ Exp[x], {х, 2, 4}, {0, 1/3, 2/3, 1, 4/3, 5/3, 2}]
Построим
график погрешности аппроксимации, то есть график разности функ
ии ril и Ехр
[х] — он представлен на рис. 11.22.
Нетрудно
заметить, что если в центральной части области аппроксимации
погрешность мала
(менее 5-10-
7
), то у правого края она резко возрастает.
Представленная
функция может использоваться и в иной форме:
Rationallnterpolation[f,{х,
m, k},{x, xmin, xmax}]
Рис.
11.22.
График погрешности рациональной аппроксимации экспоненциальной
функции
В данном
случае выбор абсцисс осуществляется автоматически в интервале от xmin до
mах.
В отличие от первого случая, когда абсциссы могли быть расположены неравномерно,
в данном случае расположение их будет равномерным. Приведем пример аппроксимации
функции синуса в интервале от n до n:
ri2 = RationalInterpolation[Sin[x],{x,3,4},{x,-Pi,Pi}]
Интересно
оценить погрешность аппроксимации. Для этого достаточно построить график разности
аппроксимирующей и аппроксимируемой функций. Это построение представлено на
рис. 11.23. Любопытно, что хотя максимальная погрешность и значительна, резких
выбросов погрешности в данном случае нет.
Рис.
11.23.
График погрешности аппроксимации синусоидальной функции
При рациональной
аппроксимации можно задать опции WorkingPrecision и Bias со значениями по умолчанию
$MachinePrecision и 0 соответственно. Опция Bias обеспечивает автоматическую
расстановку узлов интерполяции. При Bias->0 обеспечивается симметрирование
выбросов погрешности, дающее наименьшее ее значение в пиках. Ниже приведен пример
интерполяции (аппроксимации) экспоненциальной функции в интервале изменения
х от 0 до 2:
ri3
= RationalInterpolation[Exp[x],{x,2,4},{x,0,2},Bias->.25]
Построение
графика погрешности (рис. 11.24) показывает, что правильным выбором центра интерполяции
можно существенно уменьшить ее погрешность. Теперь большая погрешность наблюдается
в левой части графика. Однако резкого выброса погрешности в данном случае нет.
Рис.
11.24.
Погрешность аппроксимации экспоненты при выборе опции Bias->.25
Из приведенных
примеров ясно, что рациональная аппроксимация способна дать существенное уменьшение
погрешности при некотором оптимальном расположении узлов аппроксимации и выравнивании
погрешностей по абсолютной величине в точках минимумов и максимумов кривой погрешности.
Это лежит в основе так называемой минимаксной аппроксимации. Она реализуется
следующей функцией:
-
MiniMaxApproximation[f,{x,{xmin,xmax},m,k}]
— возвращает рациональную функцию минимаксной аппроксимации f при степени
полиномов числителя и знаменателя {m, k} ив интервале изменения х от xmin
до xmax:
-
MiniMaxApproximation
[f, approx, {x, {xmin, xmax} ,m, k} ] —возвращает рациональную функцию минимаксной
аппроксимации f при степени полиномов числителя и знаменателя
{m, k} ив интервале
изменения х от xmin до xmax с возможностью выбора метода аппроксимации
approx.
Эта аппроксимация
использует итерационный алгоритм вычислений. Они начинаются с первого шага,
на котором используется функция Rational
Interpolation. Затем аппроксимация
последовательно улучшается применением алгоритма Ремеза, лежащего в основе этого
вида аппроксимации.
Функция MiniMaxApproximation
возвращает два списка — первый с координатами абсцисс, при которых наблюдается
максимальная погрешность, второй содержит рациональную функцию аппроксимации.
Ниже представлен пример аппроксимации экспоненциальной функции:
mmlist
= MiniMaxApproximation[Ехр[х], {х, {0, 2}, 2, 4}]
Выделим формулу
аппроксимации:
mmfunc
= mmlist[[2, 1]]
Теперь можно
построить график погрешности аппроксимации (рис. 11.25).
Рис.
11.25.
График погрешности при минимаксной аппроксимации экспоненциальной
функции
Следует отметить,
что малость абсолютной ошибки для ряда функций (например, тригонометрических)
может приводить к большим относительным погрешностям в точках, где функции имеют
нулевые значения. Это может привести к отказу от выполнения аппроксимации вследствие
исчерпания числа итераций (опция Maxlterations по умолчанию имеет значение 20).
Такой случай наблюдается, например, при исполнении следующей команды:
MiniMaxApproximation[Cos[x],
{х, {1, 2}, 2, 4}]
Делением функции
на (x-Pi/2) можно исключить эту ситуацию:
MiniMaxApproximation[Cos[x]/(x-Pi/2),{*,{1!,2},2,4}]
[[2,1]]
График погрешности
для этого примера представлен на рис. 11.26. Обратите внимание на то, что в
этом примере погрешность аппроксимации не превышает (б...7)-10-
10
.
В приложении
дан список функций общей рациональной интерполяции (аппроксимации) для аналитических
зависимостей, заданных параметрически. Примеры применения этого довольно редкого
вида аппроксимации можно найти в справочной базе данных системы
Mathematica.
Там же можно найти дополнительные соображения по уменьшению погрешности аппроксимации.
Рис.
11.26.
График погрешности при минимаксной аппроксимации функции косинуса
Нули функций Бесселя — BesselZeros
В подпакете
BesselZeros определены функции, дающие список аргументов функций Бесселя в их
первых п нулевых точках: BesselJZeros [mu, n],
Bessel-YZeros[mu,n], BesselJPrimeZeros[mu,n], BesselYPrimeZeros[mu,n] и др. Ввиду редкого использования функций этого класса
ограничимся парой примеров их применения:
<<NumericalMath`BesselZeros`
BesselJZeros[0,
5]
{2.40483, 5.52008,
8.65373, 11.7915, 14.9309}
BesselJYJYZeros[2,
6/5, 3, WorkingPrecision -> 20]
{15.806622444176579073,
31.46556009153683, 47.1570167108650315}
Поиск корней уравнений с интерполяцией — InterpolateRoot
Подпакет
InterpolateRoot дает средства для ускоренного и более точного поиска корней
уравнений по сравнению с соответствующими функциями ядра. Достигается это за
счет применения интерполяции функции, корни которой ищутся. Под-пакет задает
функцию InterpolateRoot [f, {х, a, b} ], которая находит корень функции f в
интервале х от а до b. Вместо функции f можно задавать уравнение
eqn. Возможны
опции AccuracyGoal->Automatic, Maxlterations->15, WorkingPrecision->$MachinePrecision
и ShowProgress->False (указаны принятые по умолчанию значения).
Примеры применения
данной функции (n — число итераций):
<<NumericalMath`
InterpolateRoot`
n = 0; FindRoot[n++;
Exp[x] == 2, {x, 0, 1},
WorkingPrecision
-> 100, AccuracyGoal -> 95]
{x->
0.693147180559945309417232121458176568075500134360255
2541206800094933936219696947156058633269964186876}
n
17
n = 0; f[x_]
:= (n++; Exp[x]-2) /; NumberQ[x]
InterpolateRoot[f[x], {x, 0, 1), WorkingPrecision -> 100,
AccuracyGoal
-> 95]; n 10
InterpolateRoot[Exp[x]
==2, {x, 0, 1},ShowProgress -> True,
WorkingPrecision
-> 40]
{0, 0.58197670686932642439}
{21, 0, -0.12246396352039524100}
{1, 0.7019353037882764014443370764853594873432}
{21, 20, 0.0130121629575404389120930392554}
{3,0.6932065772065263165289985793736618546663}
{21, 20, 0.000062480788747713548804773113708}
{6, 0.6931471932603933841618726058237307884661}
{21, 20, 1.26443483693584888038460396742xHT8}
{12, 0.693147180559945119457822446
95590259222308309027205042483074}
{40, 20, -1.89953767048152086910014102216x
10-16}
{24, 0.6931471805599453094172321214
5786257157118117337249076750141}
Реализация интервальных методов
—IntervalRoots
Иногда важно
не найти приближенное значение корня, а уточнить интервал, в котором он находится.
В подпакете IntervalRoots для этого используется ряд известных методов, реализованных
следующими функциями:
-
IntervalBisection [f
,x, int, eps] — находит корень функции f(x) путем уточнения исходного
интервала int с заданной погрешностью eps методом половинного деления;
-
IntervalSecant [f ,x,
int, eps] — находит корень функции f(x) путем уточнения исходного
интервала int с заданной погрешностью eps методом секущей;
-
IntervalNewton [ f,
x, int, eps ] — находит корень функции/(x) путем уточнения исходного интервала
int с заданной погрешностью eps методом Ньютона (касательной).
Во всех функциях
можно опциями задать максимальное число рекурсий
(Max-Recursion) и погрешность (WorkingPrecision). Примеры применения этих функций даны ниже:
<<NumericalMath`IntervalRoots`
IntervalBisection[Sin[x], x,
Interval[{2., 8.}], .1]
Interval[{3.125,
3.218750000000001}, {6.218750000000003, 6.312500000000006}]
IntervalBisection[Sin[x],
x, Interval[{2., 8.}], .01]
Interval[{3.125,
3.17188}, {6.26563, 6.3125}]
IntervalBisection[Sin[x],
x, Interval[{2., 8.}], .01, MaxRecursion -> 10]
Interval[{3.13672,
3.14258}, {6.27734, 6.2832}]
IntervalSecant[Sin[x],
x, Interval[{2., 8.}], .01]
Interval[{3.14159,
3.1416}, {6.28316, 6.28321}]
IntervalSecant[Sin[x],
x, Interval[{2., 8.}], .01]
Interval[{3.14159,
3.1416}, {6.28316, 6.28321}]
IntervalBisection[Sin[x], x,
Interval[{2,
8}], .1, WorkingPrecision -> Infinity]
Табличное численное интегрирование — Listlntegrate
Встроенная
в ядро функция NIntegrate вычисляет определенные интегралы при известной подынтегральной
функции. Однако нередко, например при экспериментах, такая функция задается
таблицей или списком значений. В подпакете List-Integrate имеются функции для
решения этой задачи — табличного интегрирования:
-
Listlntegrate [ {yl,
y2,..., yn} ,h] — возвращает численное значение интеграла для функции, заданной
списком ординат yi при заданном шаге h по х;
-
Listlntegrate [ {yl,
y2,..., yn}, h, k] — возвращает численное значение интеграла для функции,
заданной списком ординат yi при заданном шаге h по х, используя k точек каждого
подинтервала;
-
Listlntegrate [ {{xl,
yl}, {х2, у2 },..., {хп, уп}}, k] — возвращает численное значение интеграла
для функции, заданной списком координат {х.., у.}. используя k точек
для каждого подынтервала.
Примеры применения
данной функции:
<<NumericalMath`Listlntegrate`
data = Tablet
n^2, {n, 0, 7}]
{0, 1, 4, 9,
16, 25, 36, 49}
ListIntegrate[data,
1]
343/3
Listlntegrate[{{0,0},{1,1},{2,4},{5,25},{7,49}},2]
241/2
При проведении
интегрирования для данных, заданных таблично, можно использовать интерполяцию:
арр
= Listlnterpolation[data,{{0,7}}] Integrate[app[x],{x,0,7}]
343/3
Integrate[Interpolation[{{0,0},{1,1},{2,4},
{5,25}, {7,49}},
InterpolationOrder->l][x],{x,0,7}]
241/2
Численное вычисление пределов — NLimit
В подпакете
N limit определена функция
Nlimit[expr,х->х0]
для численного
вычисления пределов выражений ехрг (см. примеры ниже):
<<NumericalMath` NLimit`
NLimit[Zeta[s]
- l/(s-l), s->l]
0.577216
N[EulerGamma]
0.577216
С помощью
команды Options [NLimit] можно просмотреть опции, которые используются функцией
NLimit по умолчанию. В этом подпакете задано также вычисление бесконечных сумм
Эйлера EulerSum[f, { i, imin, Infinity} ]. Например:
EulerSum[(-l)^k/(2k
+ 1) , {k, 0, Infinity}]
0.785398
EulerSumt(-1)^k/(2k
+1), {k, 0, Infinity},
WorkingPrecision->40,
Terms->30, ExtraTerms->30]
0.78539816339744830961566084579130322540
%- N[Pi/4,
40]
-2.857249565x
10-29
Имеется также
функция вычисления производной в численном виде:
-
ND [ f, х, хО] — вычисляет
первую производную f(x) в точке х0;
-
ND[f, {x,n} ,х0] —
вычисляет п-ю производную f(X) в точке х0. Пример вычисления
производной:
ND[Exp[Sin[x]],
х, 2]
-1.03312
Options[ND]
{WorkingPrecision->
16, Scale-> 1, Terms-> 7, Method-> EulerSum]
В некоторых
случаях вычисления могут быть ошибочными. Тогда следует использовать опции —
особенно опцию выбора метода Method. Помимо метода по умолчанию
(EulerSum) можно
использовать NIntegrate (метод интегрирования по формуле Коши).
Численное вычисление остатка — N Residue
В подпакете
NResidue имеется функция вычисления остатка NResidue
[expr, {x, x0} ] в точке
х=х0:
<<NumericalMath` NResidue`
NResidue[1/z, {z, 0}]
1. + 6.35614x
10-18 I
Residue[f, {z, 1.7}]
0
NResidue[f, {z, 1.7}]
0.259067 - 1.9353xl0-17I
l/((z+.2+.5
I)(z+.2-.5 I)) /. z -> 1.7
0.259067 +
0. I
Options[NResidue]
Обратите
внимание на возможные опции для этой функции в последнем примере.
Численное разложение в ряд — NSeries
Подпакет
NSeries вводит функцию NSeries [f, {x,xO,n}], которая дает численный ряд, аппроксимирующий
функцию f(x) в окрестности точки х = х
0
, включая термы
от (х -х
0
)
-n
до (х - х
0
)
п
.
Примеры применения
данной функции:
<<NumericalMath`NSeries`
NSeries[Sin[х],
{х, -2, 2}]
Rationalize[Chop[%]]
Rationalize[Chop[NSeries[Log[x], {x, 1, 5}, Radius -> 1/8]]]
Rationalize[Chop[NSeries[Exp[x], {x, 0, 5},
WorkingPrecision
-> 40, Radius -> 1/8]]]
Rationalize[Chop[NSeries[Exp[x], {x, 0, 5}, Radius -> 4]]]
Chop[NSeries[Zeta[s], {s, 1, 3}]]
Вычисление коэффициентов формулы интегрирования Ньютона—Котесса — NewtonCotes
Функция
NIntegrate,
имеющаяся в ядре системы Mathematica, реализует метод интегрирования Гаусса—Кронрода.
Еще одним известным методом интегрирования является метод Ньютона—Котесса,
сводящий интегрирование к вычислению взвешенных ординат функции в равномерно
расположенных точках оси абсцисс. Для реализации метода используются следующие
функции:
-
NewtonCotesWeights [n,
a, b] — возвращает список весовых коэффициентов и абсцисс узловых точек
{wi, xi} для квадратуры Ньютона—Котесса на интервале от а до
b;
-
NewtonCotesError [n,
f, a, b] — возвращает погрешность формулы Ньютона—Котесса для заданной функции
f.
Примеры применения
этих функций представлены ниже:
<<NumericalMath` NewtonCotes`
NewtonCotesWeights[5,
0, 10]
NewtonCotesError[5, f, 0, 10]
NewtonCotesError[5, f, a, a+h]
NewtonCotesWeights[5,
-0, 10, QuadratureType -> Open]
NewtonCotesError[5, f, 0, 10, QuadratureType ->
Open]
Обратите
внимание на то, что приведенные формулы готовят данные для численного интегрирования
методом Ньютона—Котесса, но не выполняют самого интегрирования.
Что нового мы узнали?
В этом уроке
мы научились:
-
Пользоваться алгебраическими
функциями пакета Algebra.
-
Применять вычислительные
функции пакета Calculus.
-
Работать с функциями
дискретной математики из пакета DiscreteMath.
-
Производить геометрические
расчеты с помощью пакета Geometry.
-
Выполнять алгебраические
вычисления с помощью пакета LinearAlgebra.
-
Пользоваться расширенными
функциями теории чисел из пакета NumberTheory.
-
Осуществлять численные
расчеты с помощью пакета NumericalMath.
|
