Численные методы: cos малых углов
нужно написать программу (Fortran) вычисления косинуса произвольного угла.
При использовании стандартной функции cos(), все значения, соответствующие углу, меньшему
0,013988227117806674 равны 1. При использовании переменных двойной точности этот диапазон лишь слегка
расширяется...
Вот это и главная проблема.
Пробовал вычислять это разложением в ряд cos(x)=1-x^2/2!+x^4/4!+...., но очередное слагаемые слишком
быстро вываливается за 1E-38...
Пробовал через синус. Тут дело гораздо лучше :) sin(x)=x-x^3/3!+x^5/5!+... . За счет sin(x)~x многое выигрывается.
Но когда пытаюсь перейти к cos(x)=sqrt(1.-sin**2), точность теряется и косинус вновь 1 :(
Быть может, кто-нибудь знает как можно решить эту проблему или хотя бы ресурсы подобного характера?
Заранее спасибо!
Вот, собственно, в чем проблема:
нужно написать программу (Fortran) вычисления косинуса произвольного угла.
При использовании стандартной функции cos(), все значения, соответствующие углу, меньшему
0,013988227117806674 равны 1. При использовании переменных двойной точности этот диапазон лишь слегка
расширяется...
Вот это и главная проблема.
Пробовал вычислять это разложением в ряд cos(x)=1-x^2/2!+x^4/4!+...., но очередное слагаемые слишком
быстро вываливается за 1E-38...
Пробовал через синус. Тут дело гораздо лучше :) sin(x)=x-x^3/3!+x^5/5!+... . За счет sin(x)~x многое выигрывается.
Но когда пытаюсь перейти к cos(x)=sqrt(1.-sin**2), точность теряется и косинус вновь 1 :(
Быть может, кто-нибудь знает как можно решить эту проблему или хотя бы ресурсы подобного характера?
Заранее спасибо!
Проблема решается только повышением точности. Есть 80-тибитные числа (поддерживаются сопроцессором, в C - long double, в Fortran - понятия не имею). Можно также попробовать считать символьно (массив символов, скажем, 300).
А дело вот в чем. cos(x) при x->0 отличается от 1 на величину 2-ого порядка малости (o(x)), в то время как sin(x)=x (при x->0), т.е. равен величине первого порядка малости. Соответственно, cos малых углов отличаются от 1 на величину o(x), что соответственно на порядок меньше, чем x. Т.е при очень малом x мы можем довольно точно принять cos(x)=1, но не можем принять sin(x)=0, так как величина x - б.м. первого порядка. Погрешность данных вычислений всегда вычисляется пропорционально аргументу. При очень большом x точность снижается, при очень малом - наоборот возрастает. Следовательно, ваша 1 при cos(x->0) гораздо точнее, чем cos((200*pi)+x->0)=1. Необходимо очень сильно увеличить точность, чтобы увидеть "все эти 10 девяток", а после них - нужные цифры. Чего обычно не делается (ну, разве что в астрономических расчетах).
Можно также увеличить число выводимых знаков после запятой, чтобы вас не пугала 1 от ненулевого аргумента cos.
Надеюсь, убедил, что cos(x) при малых x (смотрите ваше число) очень мало отличается от 1 с очень большой степенью точности. Результат, по сути верен.
Однако, я не понимаю, почему проблема возникает на уровне 1E-38? Я сам, к сожалению, уже давно не работаю с Фортраном, но, если не ошибаюсь, кроме типа REAL там были ещё и REAL*8 и даже TEMPREAL (REAL*10). При работе с такими типами минимальные представимые числа имеют порядок не 38, а 4000 с чем-то.
Посмотрите внимательно на список типов, которые поддерживает Ваша версия Фортрана, и, может быть, станете немного счастливее.
Envel всё очень классно объяснил о погрешностях. Всё действительно так и есть.
Однако, я не понимаю, почему проблема возникает на уровне 1E-38? Я сам, к сожалению, уже давно не работаю с Фортраном, но, если не ошибаюсь, кроме типа REAL там были ещё и REAL*8 и даже TEMPREAL (REAL*10). При работе с такими типами минимальные представимые числа имеют порядок не 38, а 4000 с чем-то.
Посмотрите внимательно на список типов, которые поддерживает Ваша версия Фортрана, и, может быть, станете немного счастливее.
У меня Фортрана нет. 1E-38 - это наименьшее положительное число для типа real (6 байт), все числа меньше данного считаются машинным нулем.
Видимо, где-то в вычислениях применяются real (6 байт), что и есть "узкое место".
У меня Фортрана нет. 1E-38 - это наименьшее положительное число для типа real (6 байт), все числа меньше данного считаются машинным нулем.
Видимо, где-то в вычислениях применяются real (6 байт), что и есть "узкое место".
Видимо, я неясно выразился. Прошу прощения.
У меня Фортрана тоже нет, это я намёки Tavixу даю. Пусть почитает про свой Фортран и найдёт там про вычисления с повышенной точностью. В конце концов Фортран и существует для вычислений, и не может быть, чтобы стандартный REAL там был единственным вариантом. Вот Вы писали про 80-битные числа, и я подсказываю, что в моём Фортране они назывались TEMPREAL. Это была какая-то версия компилятора от Intel... древняя очень.
Проблема решается только повышением точности. Есть 80-тибитные числа (поддерживаются сопроцессором, в C - long double, в Fortran - понятия не имею). Можно также попробовать считать символьно (массив символов, скажем, 300).
А дело вот в чем. cos(x) при x->0 отличается от 1 на величину 2-ого порядка малости (o(x)), в то время как sin(x)=x (при x->0), т.е. равен величине первого порядка малости. Соответственно, cos малых углов отличаются от 1 на величину o(x), что соответственно на порядок меньше, чем x. Т.е при очень малом x мы можем довольно точно принять cos(x)=1, но не можем принять sin(x)=0, так как величина x - б.м. первого порядка. Погрешность данных вычислений всегда вычисляется пропорционально аргументу. При очень большом x точность снижается, при очень малом - наоборот возрастает. Следовательно, ваша 1 при cos(x->0) гораздо точнее, чем cos((200*pi)+x->0)=1. Необходимо очень сильно увеличить точность, чтобы увидеть "все эти 10 девяток", а после них - нужные цифры. Чего обычно не делается (ну, разве что в астрономических расчетах).
Можно также увеличить число выводимых знаков после запятой, чтобы вас не пугала 1 от ненулевого аргумента cos.
Надеюсь, убедил, что cos(x) при малых x (смотрите ваше число) очень мало отличается от 1 с очень большой степенью точности. Результат, по сути верен.
Благодарю за столь подробное объяснение!
Хотелось бы еше немного рассказать, раз уж тема зашла так далеко...
Введение переменных типа REAL*8 не решает главной проблемы - потярю точности. Оно лишь расширяет границы допустимых значений, но всегда найдется что-то, что ни REAL, ни REAL*8 посчитать не сможет, потому как нельзя представить бесконечную периодическую дробь конечным числом знаков, пусть их будет хоть 10^1024... Это повысит точность, возможно ее будет достаточно для решения задачи, но не в этом дело... Нужно было на переменных типа REAL объяснить, почему происхоит такая вещь, как обнуление результата...
sq_deep и Envel
Спасибо большое за обсуждение темы!