2013-03-26 46 views
2

這是一個家庭作業(過去的截止日期)的觀察,我們不得不使用Euler顯式方案來研究捕食者 - 獵物模型。我比較Fortran和C代碼(下面給出),並且我無法解釋爲什麼在比較Fortran代碼和C代碼時結果存在差異。類似的C和Fortran代碼的不同結果

的C代碼:

#include <stdio.h> 

//Functions for Euler Explicit 
float FR(float R,float F,float alpha) 
{ 
    return (2*R - alpha*R*F); 
} 

float FF(float R,float F,float alpha) 
{ 
    return (-F + alpha*R*F); 
} 

int main() 
{ 
    int N; 
    float alpha,initR,initF,bigR,h; 

    /* Data Setup */ 
    alpha = 0.001; 
    initR = 15; 
    initF = 22; 
    N  = 50; 
    bigR = 400; 
    h  = 0.01; 

    float R0,F0,R1,F1; 
    int i; 

    R0 = initR; 
    F0 = initF; 

    for (i = 0; i < N; i++) 
    { 
     R1 = R0 + h*FR(R0,F0,alpha); 
     F1 = F0 + h*FF(R0,F0,alpha); 
     printf("%d\t%f\t%f\n",i,R1,F1); 
     R0 = R1; 
     F0 = F1; 
    } 

    return 0; 
} 

現在的Fortran版本,初始化參數是完全相同於上述的C版。

! EULER EXPLICIT 
SUBROUTINE eulers_explicit (initR,initF,N,alpha,h) 
    IMPLICIT NONE 
    REAL :: R0,F0,R1,F1,initF,initR,alpha,h 
    INTEGER I,N 

    R0 = initR 
    F0 = initF 

    DO I=1,N 
      R1 = R0 + h*FR(R0,F0,alpha) 
      F1 = F0 + h*FF(R0,F0,alpha) 
      PRINT *,I,R1,F1 
      R0 = R1 
      F0 = F1 
    END DO 
END SUBROUTINE eulers_explicit 
! EULER EXPLICIT 

REAL FUNCTION FR(R,F,alpha) 
    REAL, INTENT(IN)  :: F,alpha 
    REAL, INTENT(INOUT)  :: R 

    R = 2*R - alpha*R*F 
FR = R 
END FUNCTION FR 

REAL FUNCTION FF(R,F,alpha) 
    REAL, INTENT(IN)  :: R,alpha 
    REAL, INTENT(INOUT)  :: F 

    F = -F + alpha*R*F 
FF = F 
END FUNCTION FF 

PROGRAM solve 
    USE usual_routines 
    IMPLICIT NONE 
    INTEGER :: N 
    REAL :: alpha,initR,initF,h 

    alpha = 0.001 
    initR = 15 
    initF = 22 
    N  = 50 
    h  = 0.01 

    PRINT *, "Eulers explicit method: " 
    CALL eulers_explicit (initR,initF,N,alpha,h) 

END PROGRAM solve 

現在的結果,從C代碼結果的快照:

0 15.296700 21.783300 
    1 15.599301 21.568800 
    2 15.907923 21.356476 
    3 16.222683 21.146309 
    4 16.543707 20.938276 
    5 16.871117 20.732357 
    6 17.205042 20.528532 
    7 17.545610 20.326778 
    8 17.892956 20.127077 
    9 18.247213 19.929407 
    10 18.608521 19.733749 

結果從Fortran代碼:

 0 29.9666996  -21.5607319  
     1 61.1852989  20.4571381  
     2 122.330109  -18.1591854  
     3 249.350449  13.8127756  
     4 500.209259  -7.04162502  
     5 1013.98022  -2.80275159E-02 
     6 2048.26880  -2.91000959E-02 
     7 4137.56299  -9.1E-02 
     8 8358.25781  -0.668782473  
     9 16889.3262  -10.6198158  
     10 34297.5977  -353.508148 

爲什麼我看到的這種分歧? 我正在使用gfortran(v4.7.2),並在我的筆記本電腦上安裝了Intel i7上的Arch Linux。

編輯:

這是修復。

REAL FUNCTION FR(R,F,alpha) 
    REAL, VALUE :: F,alpha 
    REAL, VALUE :: R 

    R = 2*R - alpha*R*F 
FR = R 
END FUNCTION FR 

請注意VALUE用於指針分離。

+1

我建議你增加一些調試點,以便追蹤代碼中出現差異的地方。 – 2013-03-26 21:51:22

+1

爲什麼「漂浮」?如果沒有很好的理由,請始終選擇'double'作爲C中的浮點值。 – pmg 2013-03-26 21:52:39

+0

請發佈[SSCCE](http://sscce.org)。我相信你可以減少這段代碼,但仍然會重現這個問題。 – djechlin 2013-03-26 21:58:00

回答

7

FORTRAN子程序FF和FR修改它們的參數,而C等效的 沒有。

3

FF中,通過F作爲參考。在FR中,通過R作爲參考。像Fortran代碼那樣更新FR

編輯:是的,@JimLewis說。

相關問題