#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#define N 24
#define WUM 20
#define WLM 75
#define WDM 80
#define WM 175
#define Kc 0.63
#define C 0.10
#define B 0.3
#define IM 0
#define SM 20
#define EX 1
#define KG 0.3
#define KI 0.41
#define CI 0.6
#define CG 0.99
#define F 537 //模型引數輸入
int i;
double P[24]={10,24.1,20.4,18.3,10.1,5.5,0.6,3.1,1.9,4.6,5,4.8,36.2,29,6,3.6,0.4,0,0.5,3.8,0,1.8,0.2,0.3};
double Eo[24]={0.1,0,0.1,0.5,0.7,0.9,0.8,0.7,0.5,0.3,0.2,0.1,0,0,0.1,0.6,0.8,1,0.9,0.8,0.7,0.5,0.3,0.1};
double Em[N],WU[N],WL[N],WD[N],EU[N],EL[N],ED[N],E[N]; //蒸散發引數
double PE[N],Wo[N],Wmm,A[N],R[N]; //產流引數
double FRo[N],FR[N],So[N],S1[N],S2[N],SMM,SMMF,SMF,AU[N],RS[N],RI[N],RG[N],S2[N]; //三水源劃分引數
double QS[N],QI[N],QG[N],QT[N],Q[N]; //匯流引數
double UH[3]={0.3,0.6,0.1}; //無因次單位線UH的3個時段
double U;
FILE *fp;
int main ()
{
fp=fopen("計算結果.txt","w+");
WU[0]=0;WL[0]=70;WD[0]=80;Wo[0]=150;
for(i=0;i<24;i++)
{
Em[i]=Eo[i]*Kc; //蒸散發計算開始
PE[i]=P[i]-Em[i];
if (WU[i]+P[i]>=Em[i]) //上層蒸發
{
EU[i]=Em[i];
EL[i]=0;
ED[i]=0;
}
else
{
if(WL[i]>=C*WLM) //上層不夠下層補
{
EU[i]=WU[i]+P[i];
EL[i]=(Em[i]-EU[i])*WL[i]/WLM;
ED[i]=0;
}
if(C*(Em[i]-EU[i])<=WL[i])
{
EU[i]=WU[i]+P[i];
EL[i]=C*(Em[i]-EU[i]);
ED[i]=0.0;
}
else //下層不夠深層補
{
EU[i]=WU[i]+P[i];
EL[i]=WL[i];
ED[i]=C*(Em[i]-EU[i])-EL[i];
}
}
E[i]=EU[i]+EL[i]+ED[i];
fprintf(fp,"時段%d:\n(1)蒸散發引數輸出\nP=%.3f\tEo=%.3f\tEm=%.3f\tPE=%.3f\tEU=%.3f\tEL=%.3f\tED=%.3f\tE=%.3f\r\n",i+1,P[i],Eo[i],Em[i],PE[i],EU[i],EL[i],ED[i],E[i]);
Wmm=WM*(1+B); //產流計算開始
A[i]=Wmm*(1-pow(1-Wo[i]/WM,1.0/(1+B)));
if(PE[i]+A[i]<Wmm) //區域流域產流
R[i]=PE[i]-WM+Wo[i]+WM*pow(1-(PE[i]+A[i])/Wmm, 1+B);
else
R[i]=PE[i]-(WM-Wo[i]); //全流域產流
if(R[i]<0)
R[i]=0;
Wo[i+1]=P[i]-E[i]-R[i]+Wo[i];
if(Wo[i+1]>WM)
Wo[i+1]=WM;
WU[i+1]=WU[i]+Wo[i+1]-Wo[i];
WL[i+1]=WL[i];
WD[i+1]=WD[i];
if(WU[i+1]>=WUM)
{
WL[i+1]=WL[i]+WU[i+1]-WUM;
WU[i+1]=WUM;
if(WL[i+1]>=WLM)
{
WD[i+1]=WD[i]+WL[i+1]-WLM;
WL[i+1]=WLM;
if(WD[i+1]>=WDM)
WD[i+1]=WDM;
}
} //產流后流域平均蓄水容量計算
fprintf(fp,"\n(2)產流引數輸出\nA=%.3f\tR=%.3f\tW=%.3f\tWU=%.3f\tWL=%.3f\tWD=%.3f\r\n",A[i],R[i],Wo[i+1],WU[i+1],WL[i+1],WD[i+1]);
FRo[0]=0.1; //三水源劃分計算開始
if(R[i]>0)
{
So[0]=20;
FR[i]=R[i]/(PE[i]);
S1[i]=So[i]*FRo[i]/FR[i]; //初始產流面積下對應的自由水深
SMM=SM*(1+EX);
SMMF=SMM*(1-pow(1-FR[i],1.0/EX));
SMF=SMMF/(1+EX);
AU[i]=SMMF*(1-pow(1-S1[i]/SMF,1.0/(1+EX)));
if(R[i]+AU[i]>=SMMF)
{
RS[i]=FR[i]*(R[i]+S1[i]-SMF);
RI[i]=SMF*KI*FR[i];
RG[i]=SMF*KG*FR[i];
S2[i]=SMF-(RI[i]+RG[i])/FR[i]; //產流面積改變后對應的自由水深
}
else
{
RS[i]=FR[i]*(R[i]-SMF+S1[i]+SMF*pow(1-(R[i]+AU[i])/SMMF,1+EX));
RI[i]=KI*FR[i]*(R[i]+S1[i]-RS[i]/FR[i]);
RG[i]=KG*FR[i]*(R[i]+S1[i]-RS[i]/FR[i]);
S2[i]=S1[i]+R[i]-(RS[i]+RI[i]+RG[i])/FR[i];
}
FRo[i+1]=FR[i];
So[i+1]=S2[i];
}
if(R[i]==0)
{
FR[i]=0;
S1[i]=0;
SMMF=0;
SMF=0;
AU[i]=0;
RS[i]=0;
RI[i]=0;
RG[i]=0;
S2[i+1]=0;
}
FRo[i+1]=FR[i];
fprintf(fp,"\n(3)三水源劃分引數輸出\nFR=%.3f\tS1=%.3f\tSMMF=%.3f\tSMF=%.3f\tAU=%.3f\tRS=%.3f\tRI=%.3f\tRG=%.3f\tS2=%.3f\r\n",FR[i],S1[i],SMMF,SMF,AU[i],RS[i],RI[i],RG[i],S2[i]);
U=F/(3.6*2); //匯流計算開始
if(i==0)
QS[i]=U*RS[i]*UH[0];
QI[0]=40*CI+(1-CI)*RI[0]*U;
QG[0]=20*CG+(1-CG)*RG[0]*U;
if(i==1)
QS[i]=U*(RS[i]*UH[0]+RS[i-1]*UH[1]);
if(i>1)
QS[i]=U*(RS[i-2]*UH[2]+RS[i-1]*UH[1]+RS[i]*UH[0]);
if(i>0)
{
QI[i]=CI*QI[i-1]+(1-CI)*RI[i]*U;
QG[i]=CG*QG[i-1]+(1-CG)*RG[i]*U;
}
QT[i]=QS[i]+QI[i]+QG[i];
if(i==0)
Q[i]=60*UH[1]+QT[i]*UH[0];
if(i==1)
Q[i]=60*UH[2]+QT[i-1]*UH[1]+QT[i]*UH[0];
if(i>1)
Q[i]=QT[i-2]*UH[2]+QT[i-1]*UH[1]+QT[i]*UH[0];
fprintf(fp,"\n(4)匯流引數輸出\nQS=%.3f\tQI=%.3f\tQG=%.3f\tQT=%.3f\tQ=%.3f\r\n\n\n",QS[i],QI[i],QG[i],QT[i],Q[i]);
}
fclose(fp);
printf("請前往計算結果.txt查看\n");
return 0;
}
uj5u.com熱心網友回復：
暈
那你就把中間計算結果都列印出來，然后看看是哪一步計算和你手算不一樣的，再去檢查相關的邏輯處理是否有問題
uj5u.com熱心網友回復：
就是比如這個值算出來應該是15.24
但是程式給出的就是15.255
uj5u.com熱心網友回復：
浮點數又誤差，計算時要考慮誤差的疊加uj5u.com熱心網友回復：
計算精度問題。樓主筆算的時候 FR1 = 0.381
程式運行時， FR[0] = 0.3813811788774169
從而看出，程式的精度更高，計算結果比手算的更精確一些。
結果是 15.255247155096674
轉載請註明出處，本文鏈接：https://www.uj5u.com/houduan/49209.html
標籤：C語言