Цикл Фортрана не работает

Fortran возвращает сообщение об ошибке, говорящее о переполнении с плавающей запятой, когда я запускаю это

program single
implicit none
integer :: n,k
real, dimension(255) :: x
n = 255
x(1) = 1/3.0
x(2) = 1/12.0
do k = 3,n
  x(k) = 2.25*x(k-1) - 0.5*x(k-2)
end do
print *,x
end program

После некоторого тестирования я думаю, что это связано с индексом переменной в цикле do, но я не совсем понимаю, в чем проблема. Может кто-нибудь мне помочь?


fortran do-loops
person Roger    schedule 31.07.2018    source источник

Ответы (1)


arrow_upward
3
arrow_downward

Распечатав значения x во время цикла, вы можете получить некоторое представление о том, что происходит. Первые несколько значений

2.08333284E-02
5.20832092E-03
1.30205788E-03
3.25469766E-04
8.12780345E-05
2.01406947E-05
4.67754580E-06
4.54130713E-07
-1.31697880E-06
-3.19026776E-06
-6.51961273E-06
-1.30739954E-05

Обратите внимание, как меняется знак, и значение продолжает монотонно возрастать оттуда. Вы могли заметить, что для отрицательных значений x определенная вами последовательность растет экспоненциально, поэтому, как только вы станете отрицательным, экспоненциальный рост быстро доставит вас и приведет к переполнению с плавающей запятой, так как x просто растет и растет (после 100 итераций его уже за 1E+21).

Однако давайте посмотрим, что произойдет, если мы увеличим точность с плавающей запятой, указав

Integer, Parameter :: wp = Selected_real_kind( 13, 70 )
real(kind=wp), dimension(255) :: x

Если мы сейчас скомпилируем и запустим ту же самую программу, мы заметим, что x(k) становится равным 0 при увеличении k, но не меняет знак в течение 255 итераций. Так что на самом деле переполнение, которое вы видите, вызвано недостаточной точностью представления с плавающей запятой, что приводит к ошибочной смене знака и последующему экспоненциальному росту.

person Kai Guther    schedule 31.07.2018
comment
Пожалуйста, не учите (kind=8) новичков, это нестандартно и может быть не переносимым. В остальном, хороший ответ. - person Ross; 01.08.2018
comment
Интересное различие, которое я обнаружил, заключается в том, что результаты 32-битного и 64-битного компилятора различаются для одного и того же выбранного типа. Не уверен, что это компилятор, который я использую (FTN95), или другое округление между 32-битными и 64-битными инструкциями. 64-битные FTN95 и gFortran дают одинаковые результаты. Я также использовал x (k) = (4,5 * x (k-1) - x (k-2))/2, чтобы увидеть, изменило ли это округление, но без изменений. Переход на 8-байтовые вещественные числа удалил результат -ve, но разница округления 32-бит/64-бит все еще очевидна. - person johncampbell; 05.08.2018