Читать онлайн «Maple 9.5/10 в математике, физике и образовании». Страница 68

5.10.2. Аппроксимации рядом Тейлора

Начнем с аппроксимации функции хорошо известным рядом Тейлора степени 8 относительно середины интервала (точки с х=2):

> s := map(evalf, taylor(f(x), х=2, 9));

s := 0.4065945998 - 0.1565945998(x-2) + 0.00209790791(х-2)2 + 0.01762626393(х-2)3 - 0.006207547150(x-2)4 + 0.00057335662(x-2)5 + 0.00024331163(x-2)6 - 0.00010010534(x-2)7 + 0.00001414211(х-2)8 + O((x-2)9)

> TaylorApprox := convert(s, polynom):

Такой ряд позволяет использовать для вычислений только арифметические действия, что само по себе здорово! Для удобства преобразуем аппроксимацию в функцию, чтобы она соответствовала форме, указанной для первоначальной функции f(х). Тогда мы сможем построить график кривой ошибок для аппроксимации полиномом Тейлора:

> TaylorApprox := unapplу(TaylorApprox, х);

TaylorApprox := x→0.7197837994 - 0.1565945998x + 0.00209790791(x-2)2 + 0.01762626393(x-2)3 - 0.006207547150(x-2)4 + 0.00057335662(x-2)5 + 0.00024331162(x-2)6 - 0.00010010534(x-2)7 + 0.0000141421(x-2)8

Кривая ошибок для аппроксимации полиномом Тейлора строится командой

> plot(f - TaylorApprox, 0..4, color=black);

и имеет вид, представленный на рис. 5.25. Эта кривая нас, прямо скажем, не слишком радует, поскольку погрешность в сотни раз превышает заданную.

Рис. 5.25. Кривая погрешности при аппроксимации рядом Тейлора

Типичное свойство аппроксимации рядом Тейлора состоит в том, что ошибка мала вблизи точки разложения и велика вдали от нее. В данном случае, самая большая ошибка имеет место в левой оконечной точке. Чтобы вычислить значение ошибки в точке х=0, что ведет к делению на нуль (см. определение для f(х)), мы должны использовать значение предела:

> maxTaylorError := abs(limit(f(x), х=0) - TaylorApprox(0));

maxTaylorError := 0.0015029608

Итак, в самом начале наших попыток мы потерпели полное фиаско, получив совершенно неприемлемое значение погрешности в сотни раз больше заданной. Но отчаиваться не стоит, ибо, как говорят, «даже у хорошей хозяйки первый блин — комом».

Теперь опробуем рациональную аппроксимацию Паде (Pade) функции f(x) степени (4,4). Приближения, по этому разложению, будут аппроксимировать функцию более точно, и потому ошибки округления в вычислениях станут более заметными. Поэтому зададим вычисления с двумя дополнительными знаками точности:

> Digits := 12:

> s := map(evalf, taylor(f(x), x=2, 9)):

> PadeApprox := pade(s, x=2, [4,4]);

PadeApprox := (0.341034792604 + 0.0327799035348x - 0.00612783638188(x-2)2 + 0.00452991113636(x-2)3 - 0.000431506338862(x-2)4)/( 0.068484906786 + 0.465757546607x+ 0.159149610837(x-2)2 + 0.0266813683828(x-2)3 + 0.00346967791444(x-2)4)

> PadeApprox := unapply(PadeApprox, x):

Кривая ошибки для интервала [0,4] строится командой

> plot(f - PadeApprox, 0..4,color=black);

и имеет вид, показанный на рис. 5.26.

Рис. 5.26. Кривая погрешности при Паде-аппроксимации степени (4.4)

Как и при аппроксимации рядом Тейлора, ошибка здесь мала вблизи точки разложения и велика вдали от нее. Мы снова видим из графика, что для указанной функции, самая большая ошибка — в левой оконечной точке. Однако, максимальная ошибка в Паде-аппроксимации уже на порядок меньше, чем при аппроксимации полиномом Тейлора:

> maxPadeError := abs(limit(f(x), x=0) - PadeApprox(0));

maxPadeError:=0.000353777322

Это успех, показывающий, что мы на верном пути. Но, пока, погрешность остается слишком большой по сравнению с заданной.

Знатоки техники аппроксимации знают, что лучшие приближения на заданном интервале могут быть получены, используя разложение в ряд Чебышева. Это связано с тем, что ортогональные полиномы Чебышева позволяют получить аппроксимацию, погрешность которой в заданном диапазоне изменения аргумента распределена более равномерно, чем в предшествующих случаях. Выбросы погрешности на краях интервала аппроксимации в этом случае исключены

Разложим функцию f(x) на [0,4] в ряд Чебышева с точностью 1*10-8. Это означает, что все члены с коэффициентами меньше, чем эта величина, будут опущены. Такая точность обеспечивается полиномом 13 степени:

> evalf(limit(f(x), х=0));

.500000000000

> fproc := proc(x) if x=0 then 0.5 else evalf(f(x)) fi end:

> ChebApprox := chebyshev(fproc, x=0..4, 1E-8);

Можно проверить для этого примера, что кривая ошибки при аппроксимации рядом Чебышева колеблется. Поскольку ряд Чебышева был оборван на члене степени 8 (как и полином ряда Тейлора), то максимальная ошибка оказалась все еще больше заданной.

Для последующих вычислений, полезно заметить, что мы можем использовать процедуру для нахождения численных значений f(x), которая будет намного эффективнее, чем прямое определение, которое требует численного интегрирования для каждого значения х. А именно, определим процедуру численной оценки, основанную на разложении в ряд Чебышева степени 13, так как максимальная ошибка при такой аппроксимации меньше, чем 10-8, и обеспечивает для нашей цели достаточную точность. Мы определим полином Чебышева Т(х) из пакета orthopoly, и затем для эффективной оценки преобразуем его в форму Горнера:

> F := hornerform(eval(subs(T=orthopoly[T],

 ChebApprox)));

F = 0.499999998610 + (0.192405358503 + (-0.163971754264 + (-0.0083861432817 + (0.0277082269676 + (-0.00593172541573 + (-0.00132728874257 + (0.000910057654178 + (-0.000180351181100 + (0.57685696534 10-5 + ( 0.448885653549 10-5 + (-0.990274556116 10-6 + (0.925433855729 10-7 - 0.347161977631 10-8x)x)x)x)x)x)x)x)x)x)x)x)x

> F := unapply(F, x):

Схема Горнера минимизирует число арифметических операций, заменяя операции возведения в степень операциями последовательного умножения.

Теперь рассмотрим еще более точную рациональную аппроксимацию Чебышева-Паде. Это такая рациональная функция r[m, n](х) с числителем степени m и знаменателем степени n такой же, как и для разложения в ряд Чебышева. Функция r[m, n](х) согласуется с разложения в ряд ряда Чебышева f(x) членом степени m+n. Мы вычислим аппроксимацию Чебышева-Паде степени (4, 4), подобную обычной Паде-аппроксимации, успешно выполненной ранее:

> ChebPadeApprox := chebpade(F, 0..4, [4,4]);

Построим кривую ошибок:

> with(orthopoly, Т):

> plot(F - ChebPadeApprox, 0..4,color=black);

Она представлена на рис. 5.27.

Рис. 5.27. Кривая ошибки при Паде-Чебышева рациональной аппроксимации

Максимальная ошибка и на этот раз имеет место в левой оконечной точке. Величина максимальной ошибки несколько меньше, чем ошибка при аппроксимации рядом Чебышева. Главное преимущество преставления в виде рациональной функции — высокая эффективность вычислений, которая может быть достигнута преобразованием в непрерывную (цепную) дробь (см. ниже). Однако полученная максимальная ошибка чуть-чуть больше заданной:

> maxChebPadeError := abs(F(0) - ChebPadeApprox(0));

maxChebPadeError := 0.1236749 10-5

Мы достигли впечатляющего успеха и остается сделать еще один шаг в направлении повышения точности аппроксимации.

Классический результат теории аппроксимации заключается в том, что минимакс как наилучшая аппроксимация рациональной функции степени (m, n) достигается, когда кривая ошибки имеет m+n+2 равных по величине колебаний. Кривая ошибки аппроксимации Чебышева-Паде имеет нужное число колебаний, но эта кривая должна быть выровнена (по амплитуде выбросов кривой ошибки) с тем, чтобы обеспечить наилучшее минимаксное приближение. Эта задача решается с помощью функции minimax:

> MinimaxApprox := minimax(F, 0..4, [4,4], 1, 'maxerror');

MinimaxApprox :=x→ (0.174933018974 + (0.0833009600964 + (-0.02019330447644 + (0.00368158710678 - 0.000157698045886x)x)x)x)/(0.349866448284 + (0.031945251383 + (0.0622933780130) + (-0.0011478847868 + 0.0033634353802x)x)x)x)

Максимальная ошибка в аппроксимации MinimaxApprox дается значением переменной maxerror. Заметим, что мы, наконец, достигли нашей цели получения аппроксимации с ошибкой меньшей, чем 1*10-6:

> maxMinimaxError := maxerror;

maxMinimaxError := 0.585028048949 10-6

Построим график погрешности для данного типа аппроксимации:

> plot(F - MinimaxApprox,0..4,color=black);

График ошибки, представленный на рис. 5.28 показывает равные по амплитуде колебания.

Рис. 5.28. График ошибки при минимаксной аппроксимации

Таким образом, мы блестяще добились успеха в снижении погрешности до требуемого и довольно жесткого уровня. Если бы мы задались целью получить только четыре или пять точных знаков аппроксимации, что в целом ряде случаев вполне приемлемо, то могли бы получить нужный результат гораздо раньше. Нам остается оптимизировать полученную аппроксимацию по минимуму арифметических операций и проверить реальный выигрыш по времени вычислений.