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

Question

Я пытаюсь найти приближение минимаксного полинома для синуса и косинуса, используя алгоритм обмена запомнить в MATLAB.Необходима точность до 23 бит, потому что я реализую функции синуса и косинуса для плавающей запятой IEEE-754.

Используя эту ссылку здесь (см. Страницы с 8 по 15), инструкцииприведены для нахождения полинома с использованием Mathematica и Maple, однако я не уверен, как экстраполировать эти методы для MATLAB.

В соответствии с таблицей 3 мне нужно использовать полином 5-го или 6-го порядка, чтобы получить ~ 23биты (после десятичной запятой) точности.

Я планирую сначала выполнить уменьшение диапазона всех входящих тэта до значения от -pi / 4 до + pi / 4, затем выполнить функцию синуса или косинуса по мере необходимости(Конечная цель - реализовать exp (i * x) = cos (x) + i * sin (x).

Возможно, я сам смогу следовать инструкциям этой статьи, но я нездесь я знаю, как использовать функцию запоминания для моих целей. Кроме того, я не понимаю, почему автор использовал уравнение (6) (на странице 9), и не понимаю, как было определено уравнение для k (на странице 11) (где происходит из 2796201?) и почему определяющая форма многочлена, которую мы хотим в итоге изменить на sin9x) = x + kx ^ 3 + x ^ 5 * P (x ^ 2).

Будетвместо этого лучше использовать функцию firpm (так как запомнить устарела)?

Спасибо, вся помощь и руководство очень важны, а также изменения вобеспечить лучший ответ на мой вопрос.

Answer 1

Я бы не стал пытаться разработать собственные аппроксимации.Проще взять копию «Компьютерные приближения», Харт и др.У хорошей университетской библиотеки должно быть это.23 бита - это около 7 десятичных цифр, поэтому просто выберите приближение, которое даст вам необходимую точность.Вы можете выбрать либо простое полиномиальное приближение, либо использовать рациональный полином, обычно немного лучше, если вы можете терпеть это деление.

Уменьшение диапазона имеет смысл, фактически я выбрал тот же диапазон (+/ -pi / 4) в моих собственных инструментах, потому что с этим выбором диапазона особенно легко работать.

Редактировать: (Пример использования приближений можно найти в Харте.)

Здесь я найду приближение для sin (x), где x лежит в интервале [0, pi / 4].Моя цель состоит в том, чтобы выбрать аппроксимацию с абсолютной точностью не менее 1.e-7 на этом интервале.Конечно, если у вас есть отрицательное значение для x, мы знаем, что sin (x) является нечетной функцией, так что это тривиально.

Замечу, что приближения в Харте имеют тенденцию иметь форму sin (альфа пи х), где х лежит в интервале [0,1].Если бы я затем выбрал приближение для альфа = 1/2, я бы получил приближение, которое действует в течение выбранного интервала.Таким образом, для аппроксимации в интервале [0, pi / 4] мы ищем альфа = 1/4.

Далее, я ищу аппроксимацию, которая указана с абсолютной точностью не менее 7 цифр или около тогои я предпочитаю использовать рациональные полиномиальные аппроксимации, так как они имеют тенденцию быть немного более эффективными.Сканируя таблицы на странице 118 (моя копия Харта 1978 года), я нахожу приближение с альфа = 1/4, которое соответствует счету: индекс 3060.

Это приближение будет иметь вид

sin(alpha*pi*x) = x*P(x^2)/Q(x^2)

Так что теперь я перехожу на страницу, которая дает коэффициенты для SIN 3060. В моем экземпляре, он попадает на страницы 199-200.Есть 5 коэффициентов, P00, P01, P02, Q00, Q01.(Обратите внимание на несколько нестандартные научные обозначения, используемые здесь.) Таким образом, P (полином числителя) содержит 3 члена, а Q, знаменатель имеет 2 члена.Записывая это, я получаю это:

sin(alpha*pi*x) = (52.81860134812 -4.644800481954*x^3 + 0.0867545069521*x^5)/ ...
    (67.250731777791 + x^2)

Давайте попробуем это сейчас в MATLAB.

x = linspace(0,pi/4,10001);
xt = x*4/pi; % transform to [0,1]
sine = @(x) (52.81860134812*x -4.644800481954*x.^3 + ...
     0.0867545069521*x.^5)./(67.250731777791 + x.^2);

max(abs(sin(x) -sine(xt)))
ans =
   1.6424e-09

plot(x,sin(x)- sine(xt),'-')

Обратите внимание на 1e-9, прикрепленный к оси Y.

Похоже, это наиболее разумный выбор приближения для греха (x) в течение этого конкретного интервала, хотя это дает около 29 бит точности вместо 23 бит, которые вы просили.Если вы хотите выбрать другой интервал сокращения диапазона, есть несколько вариантов, которые могут сохранить термин, возможно, за несколько бит, которые вам не нужны.

log2(max(abs(sin(x) -sine(xt))))
ans =
      -29.182