Владимир Дьяконов - Maple 9.5/10 в математике, физике и образовании

Рис. 5.24. График аппроксимируемой функции

Итак, вычисление f(х) по ее интегральному представлению совершенно не эффективно. Наша цель состоит в разработке процедуры вычислений, которая дала бы 6 точных цифр результата в интервале [0..4] и требовала, по возможности, наименьшего числа арифметических операций для каждого вычисления. Втайне не вредно помечтать о том, чтобы после аппроксимации время вычислений уменьшилось бы хотя бы в несколько раз. Что получится на деле, вы увидите чуть позже. А пока войдем в дебри аппроксимации.

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

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

> TaylorApprox := convert(s, polynom):

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

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

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

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

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

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

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

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

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

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

> Digits := 12:

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

> PadeApprox := pade(s, x=2, [4,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));

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

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

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

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

> 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 := unapply(F, x):

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