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

5.10.7. Эффективная оценка рациональных функций

Полиномы числителя и знаменателя в минимаксной аппроксимации уже выражены в форме Горнера (то есть в форме вложенного умножения). Оценка полиномом степени n в форме Горнера при n умножениях и n суммированиях это наиболее эффективная схема оценки для полинома в общей форме. Однако, для рациональной функции степени (m, n) мы можем делать кое-что даже лучше, чем просто представить выражения числителя и знаменателя в форме Горнера. Так, мы можем нормализовать рациональную функцию так, что полином знаменателя со старшим коэффициентом будет равным 1. Мы можем также заметить, что вычисление рациональной функции степени (m, n) в форме Горнера требует выполнения всего m+n сложений, m+n-1 умножений и 1 деления. Другими словами, общий индекс действия есть

m + n операций умножения/деления,

m + n операций сложения/вычитания.

Вычисление рациональной функции можно значительно сократить и далее, преобразуя ее в непрерывную (цепную) дробь. Действительно, рациональная функция степени (m, n) может быть вычислена, используя только

max(m, n) операций умножения/деления,

m + n операций сложения/вычитания.

Например, если m = n, тогда эта новая схема требует выполнения только половины числа действий умножения/деления по сравнению с предшествующим методом. Для рациональной функции MinimaxApprox, вычисление в форме, выраженной выше, сводится к 9 действиям умножения/деления и 8 действиям сложения/вычитания. Число операций умножения/деления можно сократить до 8, нормализуя знаменатель к форме monic. Мы можем теперь вычислить непрерывную (цепную) дробь для той же самой рациональной функции. Вычисление по этой схеме, как это можно видеть из вывода Maple, сводятся только 4 действиям деления и 8 действиям сложения/вычитания:

> MinimaxApprox := confracform(MinimaxApprox):

> lprint(MinimaxApprox(x));

-.468860043555e-1 + 1.07858988373/

(x+4.41994160718+16.1901836591/(x+4.29118998064+70.1943521765/(x-10.2912531257+4.77538954280/(x+1.23883810079))))

5.10.8. Сравнение времен вычислений

Теперь определим время, необходимое для вычисления функции f(x) в 1000 точек, используя первоначальное интегральное определение, и сравним его с временем, требующимся для схемы MinimaxApprox в виде непрерывной дроби. Сделаем это для системы Maple 8. Так как наше приближение будет давать только 6 точных цифр, мы также потребуем 6 точных цифр и от интегрального представления функции:

> Digits := 6: st := time():

> seq( evalf(f(i/250.0) ) , i = 1..1000 ):

> oldtime := time() - st;

oldtime := 4.075

В процессе вычислений с использованием представления рациональной функции в виде непрерывной дроби иногда требуется внести несколько дополнительных цифр точности для страховки. В данном случае достаточно внести две дополнительные цифры. Итак, новое время вычислений:

> Digits := 6: st := time():

> seq( MinimaxApprox(i/250.0), i = 1..1000 ):

> newtime := time() - st;

newtime := 0.342

Ускорение вычисления при аппроксимации есть:

> SpeedUp := oldtime/newtime;

SpeedUp := 11.915205

Мы видим, что процедура вычислений, основанная на MinimaxApprox, выполняется почти в 12 раз быстрее процедуры с использованием исходного интегрального определения. Это серьезный успех, полностью оправдывающий время, потерянное на предварительные эксперименты по аппроксимации и ее оптимизации!

Заметим, что этот результат относится только к конкретному ПК и может сильно меняться при прогонке этого примера на других. Так, читатель, знакомый с учебным курсом автора по системе Maple 7 [36] обнаружит, что там в этом примере результаты были иные и куда более ошеломляющие:

oldtime := 81.805 newtime := .694 SpeedUp := 117.87464

В чем дело? А дело в том, что более ранние результаты были получены в среде Maple 7 на компьютере с процессором Pentium II с частотой 400 МГц. А новые результаты получены уже на компьютере с процессором Pentium 4 с частотой 2,6 ГГц и с системой Maple 9.5.

5.10.9. Преобразование в код ФОРТРАНа или С

Один из поводов разработки эффективной аппроксимации для вычисления математической функции заключается в создании библиотек подпрограмм для популярных языков программирования высокого уровня, таких как ФОРТРАН или С. В Maple имеются функции преобразования на любой из этих языков. Например, мы можем преобразовывать формулу для минимаксной аппроксимации в код ФОРТРАНа:

> fortran (MinimaxApprox(х));

Итак, нами показано, что правильный выбор аппроксимации для сложной функции обеспечивает уменьшение времени ее вычисления более чем на один-два порядка (!) при весьма приличной точности в 6 верных знаков и при использовании для вычислений минимального числа арифметических операций. Применение при этом средств системы Maple позволяет генерировать разложения в различные ряды, быстро вычислять рациональные аппроксимации функций и выполнять преобразования в различные специальные формы, сочетая это с мощными средствами интерактивной работы и графической визуализации, в частности с построением графиков функции и кривых ошибок при разных видах аппроксимации. Все это обеспечивает идеальную среду для решения таких задач.

5.11. Интегральные преобразования функций

5.11.1. Прямое и обратное Z-преобразования

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

Обратное Z-преобразование сводится к преобразованию комплексной функции f(z) в функцию f(z).

Эти преобразования задаются следующими функциями:

ztrans(f, n, z) — прямое преобразование функции f(n) в f(z);

invztrans(f, z, n) — обратное преобразование f(z) в f(n).

Заметим, что прямое Z-преобразование базируется на соотношении ztrans(f(n),n,z)=sum(f(n)/z^n,n=0..infinity), записанном на Maple-языке. В первых версиях системы Maple Z-преобразования выполнялись средствами библиотеки и требовали вызова командой readlib(ztrans). Но в Maple 7/8 они уже были включены в ядро системы и предварительного вызова уже не требуют. В этом убеждают следующие примеры:

> a:=ztrans(n^2,n,z);

> invztrans(a,z,n);

n²

> ztrans(cos(Pi/4*t), t, z);

> invztrans(%,z,t);

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

5.11.2. Быстрое преобразование Фурье

Преобразование Фурье широко используется в математике, физике и электрорадиотехнике. Суть этого преобразования описана чуть ниже — см. раздел 5.11.4. Ввиду широких сфер применения этого преобразования в технике часто используется его особая разновидность — быстрое преобразование Фурье или FFT (Fast Fourier Transform).

В Maple на уровне ядра реализованы функции быстрого прямого FFT и обратного iFFT преобразований Фурье для числовых данных:

FFT(m, х, у)

evalhf(FFT(m, var(x), var(y)))

iFFT(m, x, y)

evalhf(iFFT(m, var(x), var(y)))

Здесь m — целое неотрицательное число, х и у — массивы с числом элементов, кратным степени 2 (например 4, 8, 16 и т.д.), представляющие действительные и мнимые части массива комплексных чисел (данных). Функции возвращают число элементов выходных массивов, а результат преобразований помещается в исходные массивы:

> х := array([1.,2.,3.,4.]): у := array([5.,6.,7.,8.]):

> FFT(2,х,y);

4

> print(х);

[10., -4., -2., 0.]

> print(y);

[26., 0., -2., -4.]

> iFFT(2,х,y);

4

> print(x);

[1.0000000, 2.0000000, 3.0000000, 4.0000000]

> print(y);

[5.0000000, 6.0000000, 7.0000000, 8.0000000]

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

Альтернативой преобразований Фурье в наши дни стали вейвлет-преобразования. Вейвлеты это новый обширный базис для приближения произвольных зависимостей вейвлетами — «короткими» волночками разной формы, способными к масштабированию и перемещению. Вейвлеты прекрасно подходят для приближения локальных особенностей различных зависимостей, в том числе нестационарных (с параметрами, меняющимися во времени). Ознакомиться с вейвлетами и средствами работы с ними в системах MATLAB, Mathematica и Mathcad можно по книге [55]. К сожалению, в Maple готовые средства вейвлет-преобразований отсутствуют и это серьезный недостаток этих систем.