2016-10-29 2 views
1

Почему эта программа fortran создает только нули? Когда я распечатываю его, я получаю везде -0.00000! Что я сделал не так? В Matlab он работает отлично. Я не вижу причин, почему он не работает, чтобы быть честным!Простая программа, которая производит только нули, ошибка?

Похоже на свою фракцию, которая испортила ее. если я устанавливаю x равным некоторому десятичному числу, он работает.

program main 


implicit none 
    integer iMax, jMax 
    double precision, dimension(:,:), allocatable :: T 

double precision x, dx,f,L2old,L2norm,y 

integer i, j,n,bc 

n=10 

allocate(T(1:n+2, 1:n+2)) 

T=0.0d0 

do i=2,n+1 
do j=2,n+1 

x=(j+1)*1/24 
y=(i+1)*1/24 


T(i,j)= -18*(x**2+y**2)**2 

Write(*,*)'T(',i,'',j,'', T(i,j) 

end do 
end do 

Write(*,*)'T(1,1)',T(1,1) 


end program main 
+2

Я не помню * любой * Fortran, но '* 1/24' выглядит очень много, как целочисленное деление. –

ответ

3
x=(j+1)*1/24 

1/24 представляет собой целое подразделение, которое округляется до 0. Вы должны быть в состоянии заставить деление с плавающей запятой, сделав по крайней мере один из операндов с плавающей точкой, например

x=(j+1)*1.0/24.0 
1

Как указал Джим Льюис, ответ на вопрос ОП был действительно используемым целым делением.

Ничего не важно, я думаю, важно отметить, что нужно позаботиться о том, как записывается фракция с плавающей запятой. Как показывает программа OP, x был типа DOUBLE PRECISION. Тогда правильный результат должен быть

x=(j+1)*1.0D0/24.0D0 

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

В следующей программе демонстрирует проблему ::

program test 
WRITE(*,'(A43)') "0.0416666666666666666666666666666666..." 
WRITE(*,'(F40.34)') 1/24 
WRITE(*,'(F40.34)') 1.0/24.0 
WRITE(*,'(F40.34)') 1.0D0/24.0 
WRITE(*,'(F40.34)') 1.0D0/24.0D0 
end program test 

, который, как выход

0.0416666666666666666666666666666666... 
0.0000000000000000000000000000000000 
0.0416666679084300994873046875000000 
0.0416666666666666643537020320309239 
0.0416666666666666643537020320309239 

Вы четко видеть различия. Первая строка - это математический правильный результат. Вторая строка - целочисленное деление, ведущее к нулю. Третья строка показывает выход в случае, если деление вычисляется как REAL, а четвертая и пятая строки - в DOUBLE PRECISION. Пожалуйста, примите во внимание, что в моем случае REAL подразумевает 32-битное число с плавающей запятой и DOUBLE PRECISION 64-битную версию. Точность и представление как REAL, так и DOUBLE PRECISION зависят от компилятора и не определены в the Standard. Требуется только то, что DOUBLE PRECISION имеет более высокую точность, чем REAL.

4.4.2.3 Реальный тип

Реальный тип имеет значения, аппроксимирующие математические действительные числа. Процессор должен предоставить два или более методов приближения, которые определяют наборы значений для данных типа real. Каждый такой метод имеет метод представления и характеризуется значением параметра типа вида KIND. Параметр вида типа метода аппроксимации возвращается внутренней функцией KIND (13.7.89).

Если ключевое слово вещественного типа используется без параметра вида типа, то реального типа с по умолчанию реального вида задается, и значение является своего рода ВИД (0.0). Спецификатор типа DOUBLE PRECISION определяет тип real с двойной точностью; значение вида KIND (0.0D0).Десятичная точность метода точной аппроксимации с двойной точностью должна быть больше, чем у реального метода по умолчанию.

Это фактически означает, что, если вы хотите, чтобы убедиться, что ваши расчеты выполнены с использованием 32-битной, 64-битной или 128bit представления с плавающей запятой, рекомендуется использовать правильные KIND значения, как определено в собственном модуле ISO_FORTRAN_ENV.

13.8.2.21 REAL32, REAL64 и REAL128

Значения этих по умолчанию целого числа скалярных именованных констант должна быть теми параметры типа рода, которые определяют вещественный тип которого хранения размер, выраженный в битах, равен 32, 64 и 128 соответственно. Если для любой из этих констант процессор поддерживает более одного типа такого размера, зависит от процессора, какое значение предоставляется. Если процессор не поддерживает какой-либо конкретный размер, то постоянная должна быть равна -2, если процессор поддерживает виды большего размера и -1 в противном случае.

Так что это приведет к следующему коду

PROGRAM main 
    USE iso_fortran_env, ONLY : DP => REAL64 
    IMPLICIT NONE 
    ... 
    REAL(DP) :: x 
    ... 
    x = (j+1)*1.0_DP/24.0_DP 
    ... 
END PROGRAM main 
Смежные вопросы