sinX

[nikkka]

пишу калькулятор на Java. не хочу использовать стандартные библиотеки - пишу свои. как вычеслить sinX?

[caxap]

Ряд

[ИСН]

CORDIC

[Xenia1996]

nikkka в сообщении #424581 писал(а):

пишу калькулятор на Java. не хочу использовать стандартные библиотеки - пишу свои. как вычеслить sinX?

$\sin x=x-\frac{x^3}{3!}+\frac{x^5}{5!}-\frac{x^7}{7!}+\frac{x^9}{9!}+\cdots = \sum_{n=0}^\infty\frac{(-1)^nx^{2n+1}}{(2n+1)!}$

По-моему, так, если не забыла.

Такой ряд подойдёт?

Если нет, напишу другой.

-- Сб мар 19, 2011 13:13:25 --


Вот пример кода нашла: http://pasadvice.narod.ru/prog/sin.htm

[VPro]

Xenia1996 в сообщении #424653 писал(а):

nikkka в сообщении #424581 писал(а):

пишу калькулятор на Java. не хочу использовать стандартные библиотеки - пишу свои. как вычеслить sinX?

$\sin x=x-\frac{x^3}{3!}+\frac{x^5}{5!}-\frac{x^7}{7!}+\frac{x^9}{9!}+\cdots = \sum_{n=0}^\infty\frac{(-1)^nx^{2n+1}}{(2n+1)!}$

По-моему, так, если не забыла.

Такой ряд подойдёт?

Если нет, напишу другой.

-- Сб мар 19, 2011 13:13:25 --


Вот пример кода нашла: http://pasadvice.narod.ru/prog/sin.htm

Пример кода годится, только нужно добавить приведение аргумента в диапазон $[0,\pi/2]$. При использование рядов и для$ \sin$ и для $\cos$ аргумент можно привести в интервал $[0,\pi/4]$. Очень полезно для улучшения сходимости.

[nikkka]

то есть для любых Х sinX равен этому ряду?

[profrotter]

Достаточно располагать методом расчёта значений синуса на интервале $x\in[0,\frac \pi 2].$ Можно использовать рекурсивную процедуру:

$$y_k=2(1-\frac {T^2} 2)y_{k-1}-y_{k-2},$$

где $T=\frac x {N-1},$ с начальными условиями $y_0=0;y_1=T.$ Значение синуса

$$sin(x)=y_{N-1}.$$

При $N=100$ средняя абсолютная погрешность имеет прядок $10^{-5}$, при $N=1000$ средняя абсолютная погрешность имеет прядок $10^{-7}$. Рекурсия находится на границе устойчивости, поэтому потребуется тщательное тестирование программы, а большими значениями $N$ злоупотреблять не следует! В тоже время должно быть обеспечено условие $T<<1$. Достоинства метода: не нужно вычислять факториалы (длина усечённого ряда Маклорена, который приводили выше, для вычисления значений синуса на первой четверти периода должна быть не менее 20). Промежуточные значения рекурсии соответствуют значениям $y_k=sin(kT)$.

Не стоит забывать также и о возможности заранее заложить в программу рассчитанные значения синуса в 100-1000 точек на первой четверти периода, а промежуточные значения определять путём интерполирования. Лучше сплайнами (интерполяционный многочлен Лагранжа при большом количестве точек может дать сильную раскачку между узлами, кроме того становится вычислительно неэффективным). Если точек очень много можно рассмотреть возможность вообще ограничиться кусочно-линейной интерполяцией.

[Xenia1996]

nikkka в сообщении #424694 писал(а):

то есть для любых Х sinX равен этому ряду?

Для каждого вещественного $x$ его синус будет равен сумме бесконечного числа слагаемых, которые я написала. Насчёт комплексных $x$ - понятия не имею. А синус бесконечности вообще провоцирует на появление логического парадокса, поэтому он и не определён.

[Vince Diesel]

В турбопаскале это было реализовано так. Все сводится к вычислению на отрезке $[0,\pi/4]$, поскольку на $[\pi/4,\pi/2]$ можно получить с помощью костинуса и корня. А на $[0,\pi/4]$ используется многочлен наилучшего приближения 12го, вроде, порядка. Этого оказывается достаточно для точности где-то $10^{-14}$. Правда, еще надо уметь найти этот многочлен

[mihailm]

Vince Diesel в сообщении #424705 писал(а):

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

Вряд ли его вообще кто нить сможет найти точно)))
только примерно

[Vince Diesel]

Ну конечно, надо знать его коэффициенты с точностью $10^{-16}-10^{-18}$, чтобы считать численно.

[caxap]

А как карманные калькуляторы его считают?

[VPro]

Vince Diesel в сообщении #424705 писал(а):

В турбопаскале это было реализовано так. Все сводится к вычислению на отрезке $[0,\pi/4]$, поскольку на $[\pi/4,\pi/2]$ можно получить с помощью костинуса и корня. А на $[0,\pi/4]$ используется многочлен наилучшего приближения 12го, вроде, порядка. Этого оказывается достаточно для точности где-то $10^{-14}$. Правда, еще надо уметь найти этот многочлен

Вычисление полинома 12-го порядка в общем случае требует 22 операции умножения и 12 операций сложения.
С другой стороны, если мы просто возьмем разложение в ряд до 15 степени, заранее вычислим обратные факториалы, то получим полином, для вычисления которого нужно 15 операций умножения и 7 операций сложения (за счет наличия только нечетных степеней). А точность при $x\in[0,\pi/4]$ будет, по крайней мере, не не хуже.

[ewert]

Vince Diesel в сообщении #424705 писал(а):

используется многочлен наилучшего приближения 12го, вроде, порядка. Этого оказывается достаточно для точности где-то $10^{-14}$. Правда, еще надо уметь найти этот многочлен

Найти наилучший многочлен не просто, а очень просто. Вот код на Матлабе построения многочлена 12-й степени для отрезка $0;\frac{\pi}{2}]):

Используется синтаксис Matlab M

a=0;    b=pi/2;
n=11;   k=0:n;
% Стандартные чебышёвские узлы:
z=cos((pi+2*pi*k)/(2*n+2));
% Перемасштабируем на нужный отрезок:
x=(a+b)/2 - z*(b-a)/2;

% Корректируем вручную первый узел:
x(1)=x(1) + (x(2)-x(1))*0.06;
%x(3)=x(3) + (x(3)-x(2))*0.03;

% Находим коэффициенты многочлена:
p=polyfit(x, sin(x)./x, n);

% Отладочный вывод на экран:
t=a+1e-8:(b-a)/1000:b;
r=sin(t)./t - polyval(p, t);
plot(t,r, x,x*0,'o'),  grid,  pause

% Это пока был sin(x)/x. Умножаем на икс:
p=[p,0]'

% Отладочный вывод на экран:
t=a+1e-8:(b-a)/1000:b;
r=sin(t) - polyval(p, t);
plot(t,r, x,x*0,'o'),  grid

Результат:

$P_{12}(x)=1.369970812617591\cdot10^{-9}\,x^{12}-3.041501739019810\cdot10^{-8}\,x^{11}+1.231369474597185\cdot10^{-8}\,x^{10}+2.737603227735816\cdot10^{-6}\,x^{9}+1.774771739477244\cdot10^{-8}\,x^{8}-1.984243441621931\cdot10^{-4}\,x^{7}+5.051048313043363\cdot10^{-9}\,x^{6}+8.333331938139212\cdot10^{-3}\,x^{5}+2.289140935408664\cdot10^{-10}\,x^{4}-1.666666666863060\cdot10^{-1}\,x^{3}+6.781558769758184\cdot10^{-13}\,x^{2}+9.999999999999953\cdot10^{-1}\,x^{1}$

Максимальная относительная погрешность -- где-то $4\cdot10^{-15}$. Поскольку многочлен всё-таки не совсем наилучший, этот результат можно улучшить где-то процентов на 20..30, но не более того, и для этого уже придётся потрудиться; а стоит ли?... Другое дело, что что этот способ получения многочлена -- неграмотный. Слагаемые очень сильно различаются по масштабу, что приводит к увеличению погрешностей округления и при построении многочлена, и при его использовании. По-хорошему надо было бы строить интерполяционную систему уравнений вручную, беря в качестве базисных многочлены Чебышёва и вычисляя значения полученного многочлена тоже не в явном виде, а в виде разложения по многочленам Чебышёва. Но это уже требует некоторого действительно программирования, хотя и несложного. Впрочем, результат окажется не шибко лучше -- для данного набора параметров. Вот если искать, скажем, многочлен 15-й степени, то Матлаб собственными средствами уже не справится, и придётся действительно программировать самостоятельно.

На всякий случай -- распечатка результата (вдруг я чего напортил при редактировании):

p =

    1.369970812617591e-009
   -3.041501739019810e-008
    1.231369474597185e-008
    2.737603227735816e-006
    1.774771739477244e-008
   -1.984243441621931e-004
    5.051048313043363e-009
    8.333331938139212e-003
    2.289140935408664e-010
   -1.666666666863060e-001
    6.781558769758184e-013
    9.999999999999953e-001
                         0