Файл: Mathematica для математиков. Часть Реализация основных понятий математического анализа.pdf

131
????????????????????????????????????????[???????? + ???????? ≤ ???????? && ???? + ???????? ≥ ???? && ???? ≥ ???? && ???? ≥ ????,
{????, −????, ????}, {????, −????, ????}]
В некоторых случаях функции NMaximize и NMinimize могут заменять функции Solve и NSolve. В следующем примере функция NMaximize находит единственный вещественный корень уравнения
????????????????????????????????????[{????, ????
????
+ ???? + ???? == ????}, {????}]
{1.,{x
-0.68232781}}
Сравните
????????????????????????[????
????
+ ???? + ???? == ????, ????]
{{???? → −0.68232781},
???? → 0.3411639 − 1.1615414ⅈ ,
{???? → 0.3411639 + 1.1615414ⅈ}}
Для функций NMaximize и NMinimize характерно использование опций, которые уместно использовать в численных алгоритмах.
Опция WorkingPrecision задает количество значащих цифр, используемых во внутренних вычислениях; установка WorkingPrecision->
MachinePrecision (по умолчанию) приводит к вычислениям с процессорной точностью.
Опция AccuracyGoal->n задает количество значащих цифр абсолютной погрешности, которая считается равной
n

10 . Опция
PrecisionGoal->m определяет «относительную погрешность»
m
x


10
вычислений и задается количеством значащих цифр. Установкой по умолчанию для этих опций является WorkingPrecision/2. Например, команда
????????????????????????????????????[????????????[????], ????, ???????????????????????????????????????????????? → ????????, ???????????????????????????????????????????????????? → ????]
{-1.,{x
-1.5707963}} определяет следующий критерий остановки вычислений (погрешность)


x
x
x
8 12
*
10
,
10
max




и
 
12 10



x
f
. Первое неравенство относится к координатам точки экстремума, второе – к значению целевой функции.
Опцией EvaluationMonitor задается набор команд, которые вычисляются каждый раз, когда вычисляется значение выражения, стоящего в функции NMinimize.
Пример. Найдем максимальное значение функции
 


2 2
sin
,
y
x
x
y
x
f


в области
2
/
2
/





x
,
2
/
2
/





y
. Вначале построим график функции.
????[????_, ????_]: = ????????????????[????
????
+ ????
????
]
????????????????????????[????[????, ????], {????, −????/????, ????/????}, {????, −????/????, ????/????}, ???????????????????????????????????? → ????????????]

132
Мы видим, что есть экстремумы в угловых точках области и локальные экстремумы внутри. Найдем их.
????????????????????????????????[{????[????, ????], −????/???? ≤ ???? <= ????/2, −????/???? ≤ ???? <= ????/2}, {????, ????}]
{−
1 2
????Sin[
1 4
(6 − ????)???? +
????
2 4
], {???? → −
????
2
, ???? → −
1 2
(6 − ????)????}}
????[%]
{1.5707963,{x
-1.5707963,y-1.4983284}}
Функция Maximize нашла максимум в области (один из двух имеющихся).
Теперь посмотрим, что даст функция NMaximize.
????????????????????????????????????[{????[????, ????], −????/???? ≤ ???? <= ????/2, −????/???? ≤ ???? <= ????/2}, {????, ????}]
{1.3076194,{x
1.3552111,y3.010915·10
-9
}}
Функция NMaximize нашла локальный максимум внутри области.
Изобразим последовательность «пробных» точек, которые использовала функция NMaximize при приближении к локальному максимуму. Вначале запомним координаты «пробных» точек, которые расположены недалеко от экстремума. Для этого можно использовать опцию EvaluationMonitor.
{????????????, ????????????} = ????????????????[????????????????????????????????????[{????[????, ????], −????/???? ≤ ???? ≤ ????/????, −????/???? ≤ ???? ≤ ????/????},
{????, ????}, ???????????????????????????????????????????????????????????????????? ⧴ ????????????[{????, ????}]]];
Напомним, что значение опции EvaluationMonitor является выражение, которое вычисляется каждый раз, когда вычисляется значение целевой функции. Функция Sow «выбрасывает в окружающее пространство» результат вычисления своего аргумента (в нашем случае пару {x,y}, содержащую координаты точки). Эти пары из «окружающего пространства подбирает» ближайшая объемлющая функция Reap. Первым элементом она возвращает результат вычисления основного выражения (в нашем случае результат помещается в переменную sol). Вторым элементом, который она возвращает, является список «собранных» ею пар {x,y} координат «пробных» точек.
Можете посмотреть, чем являются переменные sol и pts.
????????????
{1.3076194,{x
1.3552111,y3.010915·10
-9
}}
????????????
{{{1.4894175,1.0840092},{-1.0709764,-0.44463881},...}
Мы не стали здесь приводить длинный список координат «пробных» точек.
Теперь, имея координаты «пробных» точек, выделим и изобразим только те, которые находятся рядом с экстремумом.
????????????????????????????????????????????[????[????, ????], {????, −????/????, ????/????}, {????, −????/????, ????/????},
???????????????????????????????????????????????????? → ???????????? ????. ????# & , ???????????????????????????????? → ????????,
???????????????????????? → {????????????????????, ????????????????????????????????????[????. ????????], ????????????[????????????????????,
????????????????????[????????????????????[????????????], ????_/; ????????????[????@@???? − ????????????????????[????????????]] ≤ ????. ????]]}]

133
Здесь выражение First[pts] имееет вид
????????????????????[????????????]
{{1.4894175,1.0840092},{-1.0709764,-0.44463881},...}
Функция Cases[{e
1
,e
2
,...},шаблон] возвращает список элементов e i
, которые соответствуют шаблону. В нашем случае элементами являются пары чисел {x
i
,y
i
}. В результате шаблон z_ представляет список таких пар (т.е. он имеет вид List[x
i
,y
i
]). Условная операция /; (слэш – точка с запятой) выделяет из всех пар {x
i
,y
i
} только те, которые удовлетворяют условию, записанному после операции. В нашем случае выделяются только те точки/списки {x
i
,y
i
}, в которых значение целевой функции f[x
i
,y
i
] находится близко к найденному экстремальному значению First[sol]. Само значение f[x
i
,y
i
] получается использованием операции @@ (функции Apply), которая в нашем случае действует следующим образом: f@@z=Apply[f,{xi,yi}]= f[xi,yi].
Напомним, что функция Apply[f,expr] заменяет именем f заголовок выражения expr (в нашем случае List). С другими функциями, использованными в коде, мы знакомились ранее.
□
В функциях NMaximize и NMinimize можно использовать несколько разных алгоритмов численного поиска экстремума, которые можно задать с помощью опции Method. Иногда выбор одного из них, отличного от выбора по умолчанию, может дать лучший результат. Мы не ставим перед собой задачу изложения этих методов. С ними вы можете познакомиться по справочной системе.
У функций NMaximize и NMinimize есть опция MaxIterations , само название которой говорит о ее назначении, и опция StepMonitor, с которой мы познакомимся позже.
Для поиска локального минимума/максимума функции
 
x
f
, в окрестности значения
0
x
x

, используются функции FindMinimum и FindMaximum.
Они имеют идентичные форматы вызова, поэтому для краткости мы будем описывать форматы только одной из них. Результат возвращается в виде таком же, в каком возвращают результат функции NMaximize и
NMinimize.
В формате FindMinimum[f[x],{x,x
0
}] численно ищется локальный минимум выражения f[x]; начальное значение задается точкой x
0
. В формате
FindMinimum[f[x],x] численно ищется локальный минимум выражения f[x], начиная поиск из случайно выбираемой точки.

134
В формате FindMinimum[f[x,y,...],{{x,x
0
},{y,y
0
},...}] численно ищется локальный минимум функции нескольких переменных в окрестности точки {x
0
,y
0
,...}. В формате
FindMinimum[{f[x,y,...],ограничения},{{x,x
0
},{y,y
0
},...}] численно ищется локальный минимум с учетом ограничений. При этом точку {x
0
,y
0
,...} можно не задавать – она будет автоматически выбираться в области, определяемой ограничениями.
В следующем коде мы вызываем функцию FindMinimum три раза с разными начальными приближениями. При этом в третьем вызове мы задаем еще интервал изменения независимой переменной. Положение минимума, каждый раз разное, показано на графиках.
???? ????_ = ???? ????????????[????];
???????? = ????????????????????????????????????????????[????[????], {????, ????}]
???????? = ????????????????????????????????????????????[????[????], {????, ????????}]
???????? = ????????????????????????????????????????????[{????[????], ???? ≤ ???? ≤ ????????}, {????, ????}]
????????????????????[????????_] ≔ ????????????????[????[????], {????, ????, ????????}, ???????????????????????? → {???????? = ????/. ????????[[????]];
{????????????, ????????????????????????????????????[????. ????????], ????????????????????[{????????, ????[????????]}]}}];
{????????????????????[????????], ????????????????????[????????], ????????????????????[????????]}
{-4.8144699,{x
4.9131804}}
{-11.040708,{x
11.085538}}
{-1.6764929,{x
6.}}
Как видим, выбор начальной точки влияет на положение точки минимума. В третьем вызове ближайший к начальной точке минимум оказался на границе области.
Аналогично для функции многих переменных, стартуя из разных точек, получаем разные локальные экстремумы. Рассмотрим функцию график которой показан на следующем рисунке.
????????????????????????[???????????? ???? ????????????[????????], {????, −????, ????}, {????, −????, ????}]
Вызовем функцию FindMinimum с разными начальными точками.
????????????????????????????????????????????[????????????[????]????????????[????????], {{????, ????}, {????, ????}}]
{-1.,{x
1.5707963,y2.3561945}}
????????????????????????????????????????????[????????????[????]????????????[????????], {{????, ????}, {????, ????}}]
{-1.,{x
-1.5707963,y0.78539816}}
При задании области ограничений переменных начальную точку задавать не обязательно

135
????????????????????????????????????????????[{????????????[????]????????????[????????], ????
????
+ ????
????
≤ ????}, {????, ????}]
{-0.66871014,{x
-0.79512351,y0.60644782}}
На следующем рисунке показана целевая функция в области единичного круга, который использовался при поиске экстремума.
????????????????????????[????????????[????]????????????[????????], {????, −????. ????, ????. ????}, {????, −????. ????, ????. ????},
???????????????????????????????????????????????????????? → ????????????????????????????????[{????, ????, ????}, ????
????
+ ????
????
≤ ????]]
Ограничения могут содержать условие «целочисленности» некоторых (или всех) переменных. Зададим следующую область
????????????????????????????????????????[???????? + ???????? ≤ ???????? && ???? + ???????? ≥ ???? && ???? ≤ ???????? && ???? ≥ ????,
{????, −????, ????}, {????, −????. ????, ????}, ???????????????????????????????????????????? → ????????????????????????????????????]
Теперь найдем минимум функции x+y в этой области, задавая условие
«целочисленности» переменной y, и не задавая его.
????????????????????????????????????????????[{???? + ????, ???????? + ???????? ≤ ???? && ???? + ???????? ≥ ???? && ???? ≤ ???????????? &&
???? ≥ ???? && ???? ∈ ????????????????????????????????}, {????, ????}]
{3.3,{x
0.3,y3}}
????????????????????????????????????????????[{???? + ????, ???????? + ???????? ≤ ???? && ???? + ???????? ≥ ???? && ???? ≤ ???????? ???? &&
???? ≥ ????}, {????, ????}]
{3.1428571,{x
0.28571429,y2.8571429}}
Естественно, что результат оказался разный.
Для функций
FindMinimum и
FindMaximum характерно использование опций, которые уместно использовать в численных алгоритмах.
Опции
WorkingPrecision,
AccuracyGoal,
PrecisionGoal используются точно также, как было описано для функций NMaximize и
NMinimize
. Опция Method имеет другие значения, поскольку алгоритмы поиска локального экстремума отличаются от алгоритмов поиска глобального.
Новым является добавление опции Gradient. Она задает выражение для вычисления градиента минимизируемой / максимизируемой функции.
????[????_, ????_]: = ????????????[????
????
− ????
????
] + ????????????[????
????
+ ????
????
];
????????????????????????[????[????, ????], {????, −????, ????}, {????, −????, ????}]

136
????????????????????????????????????????????[????[????, ????], {????, ????}, ???????????????????????????????? ⧴
{????????????????????[????
????
+ ????
????
] − ????????????????????[????
????
− ????
????
], ????????????????????[????
????
+ ????
????
] + ????????
????
????????????[????
????
− ????
????
]}]
{-2.,{x
1.3446341,y1.7042148}}
Опция StepMonitor позволяет контролировать шаги вычислений.
????[????_] = ????
????
+ ????/ ???? − ????;
????????????????????????????????????????????[????[????], {????, ????. ????}, ???????????????????????????????????????????? ⧴ ????????????????????["???? = ", ????]]
x = 0.57904634 x = 0.573647 x = 0.57435347 x = 0.57434918
{0.64938489,{x
0.57434918}}
Можно, например, контролировать количество итераций
???? = ????; {????????????????????????????????????????????[????[????], {????, ????. ????}, ???????????????????????????????????????????? ⧴ ???? + +], ????}
{{0.64938489,{x
0.57434918}},4}
Используя эту опцию, можно изобразить промежуточные точки на каждом шаге вычислений.
{????????????, ????????????} = ????????????????[????????????????????????????????????????????[????[????], {????, ????},
???????????????????????????????????????????? ⧴ ????????????[{????, ????[????]}]]]
????????????????[????[????], {????, ????. ????, ????}, ???????????????????????? → {????????????, ????????????????????????????????????[????. ????????],
????????????[????????????????????, ????????????[[????]]]}, ???????????????????????????????????? → {????, ????}]
У опций EvaluationMonitor и StepMonitor есть различие в моменте вызова. Первая вызывается после каждого вычисления значения входной функции, вторая – после каждого шага итераций.

137
На рисунке в середине показан контурный график исследуемой функции и ломаная, вершины которой представляют точки последовательного приближения. Желтая точка показывает положение максимума, который определила функция FindMaximum.
????????????????????????????????????????????[????[????, ????], {????, −????, ????}, {????, −????, ????}, ???????????????????????????????????????????????????? → (????????????[????. ????#]&),
???????????????????????????????? → ????????, ???????????????????????????????????????????????????????? → ???????????????????????????????? ????, ????, ???? , ????
????
+ ????
????
≤ ???? ,
???????????????????????? → {????????????, ????????????????[????????????], ????????????????????????????????????[????. ????????], ????????????????????[????????????], ????????????????????????,
????????????????????????????????????[????. ????????], ????????????????????[{????????, ????????}]}](* пред. рисунок в центре *)
Для сравнения изобразим ломаную, по которой функция NMaximize приближается к максимуму.
{????????????????, ????????????????} = ????????????????[????????????????????????????????????[{????[????, ????], ????
????
+ ????
????
≤ ????}, {????, ????},
???????????????????????????????????????????? ⧴ ????????????[{????, ????}]]];
{????????????, ????????????} = {????, ????}/. ????????????????[[????]];
????????????????????????????????????????????[????[????, ????], {????, −????, ????}, {????, −????, ????},
???????????????????????????????????????????????????? → ???????????? ????. ????# & , ???????????????????????????????? → ????????,
???????????????????????????????????????????????????????? → ???????????????????????????????? ????, ????, ???? , ????
????
+ ????
????
≤ ???? ,
???????????????????????? → {????????????, ????????????????[????????????????], ????????????????????????????????????[????. ????????], ????????????????????[????????????????],
????????????????????????, ????????????????????????????????????[????. ????????], ????????????????????[{????????????, ????????????}]}]
Полученный контурный график представлен на предыдущем рисунке справа.
Напомним, что начальное приближение в области ограничений (если оно не задано) функция NMaximize выбирает самостоятельно.
□
3.7
Приближение функций.
Когда функция задана в дискретном наборе точек и требуется определить ее значение в промежуточных точках, то прибегают к интерполяции или аппроксимации. При этом неизвестная истинная функция заменяется некоторой другой, часто кусочной, имеющей достаточно простые выражения на кусках.
Если приближающая функция принимает заданные значения в узлах, то говорят об интерполяции, если нет, то говорят об аппроксимации.

138
3.7.1
Полиномиальная интерполяция
Функция InterpolatingPolynomial создает полином (выражение, степенной многочлен), который проходит через узловые точки. В одномерном случае по n точкам строится полином степени n-1.
При вызове функции InterpolatingPolynomial ей передается список координат узлов и идентификатор переменной полинома, например, x, поскольку возвращаться будет выражение, а не функция. При этом интерполяционный полином возвращается в форме Горнера. Например,
????????????????????????????????????????????????????????????????????????????????????????????[{????, ????, ????????, ????????, ????????????}, ????]
1 + (−1 + ????)(7 + (−2 + ????)(3 + ????))
Здесь список {1,8,27,64,125} представляет ординаты узловых точек; их абсциссами полагаются натуральные числа {1,2,3,...}. Результирующий полином всегда можно преобразовать к обычному виду.
????????????????????????[%]
????
3
Существует бесчисленное множество полиномов, проходящих через заданные точки. Функция InterpolatingPolynomial пытается возвратить полином наименьшей степени.
В формате InterpolatingPolynomial[{{x
1
,y
1
},{x
2
,y
2
},...},x] задаются обе координаты узловых точек.
В следующем коде мы выполняем полиномиальную интерполяцию по множеству точек, приближенно расположенных на параболе
2
x
y

; смещение с параболы задается небольшим случайным значением.
????????????????????????????????????????[ ];
???????????????? = ???????????????????? ????, ????
????
+ ???????????????????????????????????????? −????. ????, ????. ???? , ????, −????, ????
????[????_] = ????????????????????????????[????????????????????????????????????????????????????????????????????????????????????????????[????????????????, ????], ????]
????????????????[????[????], {????, −????, ????}, ???????????????????????????????????? → {????????????????????, ????????????????????????????????????[????. ????????]},
????????????????????????−> {????????????, ????????????????????????????????????[????. ????????], ????????????????????/@????????????????}]
-0.29820805+0.33587531 x+1.9489638 x
2
-0.18505705 x
3
-
0.42240753 x
4
+0.029138328 x
5
+0.061215827 x
6
-0.0017099349 x
7
-
0.0034036856 x
8
+0.000032814108 x
9
+0.000062829842 x
10
Здесь функция SeedRandom[ ] инициализирует генератор случайных чисел текущим временем.
В следующем пункте мы рассмотрим кусочно – полиномиальную интерполяцию.
Забегая вперед, выполним сравнение кусочно – полиномиальной интерполяции, выполняемой функцией Interpolation, и полиномиальной, выполняемую функцией InterpolatingPolynomial.
???????????????? = ????, ???? , ????, ???? , ????, ???? , ????, ???? , ????, ???? ;
???? = ???????????????????????????????????????????????????? ???????????????? ;

139
????[????_] = ????????????????????????????[????????????????????????????????????????????????????????????????????????????????????????????[????????????????, ????], ????]
????????????????[????????????????[{????[????], ????[????]}, {????, ????, ????}, ???????????????????????????????????? → {????????????????????, ????????????????}],
???????????????????????????????? ????????????????, ???????????????????????????????????? → ????????????, ???????????????????????????????????? ????. ???????? ,
???????????????????????????????????? → {{−????, ????}, {−????, ????}}, ???????????????????????????????????????? → {????, ????}]
1 −
11????
30
+
683????
2 120
−
79????
3 30
+
37????
4 120
3.7.2
Кусочно – полиномиальная интерполяция
Для представления результатов интерполяции в Mathematica используется специальный объект InterpolatingFunction. Это «чистая» функция, которую можно использовать как обыкновенную функцию: вычислять значение в точке, дифференцировать и т.д. Единственное ограничение состоит в том, что она создается на определенном диапазоне изменения аргументов. При попытке подстановки аргумента, выходящего за границы этого диапазона, генерируется сообщение.
InterpolatingFunction выполняет кусочно – полиномиальную интерполяцию, порядок гладкости которой можно задать при создании функции. Она обычно создается другими функциями системы, выполняющими приближенные вычисления. В частности при численном решении ДУ с помощью функции NDSolve создается приближенное решение в виде такой функции.
Алгебраическую интерполяцию в системе выполняет функция
Interpolation
, которая принимает набор точек и создает объект
InterpolatingFunction.
???????????????? = {????, ????, ????, ????, ????, ????};
???? = ????????????????????????????????????????????????????[????????????????]
???? ????. ????
????????????????[????????????????[????[????], {????, ????, ????}, ???????????????????????????????????? → {????????????????????, ????????????????????????????????????[????. ????????????]}],
???????????????????????????????? ????????????????, ???????????????????????????????????? → ????????????, ???????????????????????????????????? ????. ???????? ,
???????????????????????????????????? → {{????, ????}, {−????, ????}}, ???????????????????????????????????????? → {????, ????}]
InterpolatingFunction[{{1,6}}, " <> "]
2.25
Здесь мы использовали формат вызова Interpolation[{y
1
,y
2
,...}], который полагает, что абсциссы точек равны {1,2,...}, а ординаты –
{y
1
,y
2
,...}.