這是一個家庭作業(過去的截止日期)的觀察,我們不得不使用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
用於指針分離。
我建議你增加一些調試點,以便追蹤代碼中出現差異的地方。 – 2013-03-26 21:51:22
爲什麼「漂浮」?如果沒有很好的理由,請始終選擇'double'作爲C中的浮點值。 – pmg 2013-03-26 21:52:39
請發佈[SSCCE](http://sscce.org)。我相信你可以減少這段代碼,但仍然會重現這個問題。 – djechlin 2013-03-26 21:58:00