# 一维对流扩散问题TDMA算法

/*

upwind –left=(左端点处的值) –right=(右端点处的值) –pdelta=(网格Pe数)
–scheme=(格式，可选center/upwind/mixed/exponent/power/precise) –output=(输出文

upwind –left=100 –right=200 –pdelta=5 –scheme=upwind –output=Pe5Upwind.dat

//*/
#include <stdio.h>
#include <math.h>
#include <string.h>
#include <stdlib.h>

double max(double number1,double number2);
void parseParameter(int *argc, char **argv, double *PHI_0, double *PHI_L, double *P_Delta, char *scheme, char *outputfilename);

int main(int argc, char **argv)
{

const int GRIDAMOUNT=20;//网格数目为20，由于两端点已知，实际未知数个数为19个。
double PHI_0=100;

double PHI_L=200;
double P_Delta=1; //网格Pe数，取为1,5,10等。
char scheme[10]="upwind";

int PRECISE=0;
char outputfilename[100]="output.dat";

double a_E[GRIDAMOUNT-1];
double a_W[GRIDAMOUNT-1];

double a_P[GRIDAMOUNT-1];
double A_P_Delta[GRIDAMOUNT-1];

double P_TDMA[GRIDAMOUNT-1];
double Q_TDMA[GRIDAMOUNT-1];

double PHI[GRIDAMOUNT-1];

int i;

FILE *dataoutput;

parseParameter(&argc, argv, &PHI_0, &PHI_L, &P_Delta, scheme, outputfilename);

if(!strcmp(scheme,"upwind"))
{
for (i=0;i<GRIDAMOUNT-1;i++){A_P_Delta[i]=1;}

}
else if(!strcmp(scheme,"center"))

{
for (i=0;i<GRIDAMOUNT-1;i++){A_P_Delta[i]=1-0.5*fabs(P_Delta);}

}
else if(!strcmp(scheme,"mixed"))

{
for (i=0;i<GRIDAMOUNT-1;i++){A_P_Delta[i]=max(0,1-0.5*fabs(P_Delta));}

}
else if(!strcmp(scheme,"exponent"))

{
for (i=0;i<GRIDAMOUNT-1;i++){A_P_Delta[i]=fabs(P_Delta)/(exp(fabs(P_Delta))-1);}

}
else if(!strcmp(scheme,"power"))

{
for (i=0;i<GRIDAMOUNT-1;i++){A_P_Delta[i]=max(0,pow((1-0.1*fabs(P_Delta)),5));}

}
else if(!strcmp(scheme,"precise"))

{
PRECISE=1;
}
else
{

fprintf(stderr,"\nUnknown scheme %s!\nAccepted schemes are: upwind|center|mixed|exponent|power|precise\n",scheme);
exit(2);

}
//      A_P_Delta[i]=1-0.5*fabs(P_Delta);

//      A_P_Delta[i]=1;
//      A_P_Delta[i]=max(0,1-0.5*fabs(P_Delta));

//      A_P_Delta[i]=fabs(P_Delta)/(exp(fabs(P_Delta))-1);

//      A_P_Delta[i]=max(0,pow((1-0.1*fabs(P_Delta)),5));

/*
这里A(|P_Delta|)的取法可见陶书第一版P201，各种格式的取法如下：
中心差分：1-0.5*fabs(P_{\Delta})
迎风差分：1
混合格式：max(0,1-0.5*fabs(P_{\Delta}))
指数格式：fabs(P_{\Delta})/(exp(fabs(P_{\Delta}))-1)
乘方格式：max(0,(1-0.1*fabs(P_{\Delta}))^5)
注意网格P_Delta和整体的P_e的不同！

参考精确解：(exp(P_Delta*(i+1))-1)/(exp(P_Delta*GRIDAMOUNT)-1)*(PHI_L-PHI_0)+PHI_0
*/

/* 下面开始对诸格式均一致的计算过程 */

/* a_W,a_E,a_P赋值 */
for (i=0;i<GRIDAMOUNT-1;i++)

{
a_E[i]=A_P_Delta[i]+max(0,-P_Delta);//此处P_Delta在东西界面上相同，所以不区分。
a_W[i]=A_P_Delta[i]+max(P_Delta,0);

a_P[i]=a_E[i]+a_W[i];//后边的+Fe-Fw项由常物性假定而略去
}

/* TDMA算法求解这一方程组 */
P_TDMA[0]=a_E[0]/a_P[0];

Q_TDMA[0]=a_W[0] * PHI_0 / a_P[0];

for (i=1;i<GRIDAMOUNT-2;i++)

{
P_TDMA[i]=a_E[i]/(a_P[i]-a_W[i]*P_TDMA[i-1]);

Q_TDMA[i]=(a_W[i]*Q_TDMA[i-1])/(a_P[i]-a_W[i]*P_TDMA[i-1]);

}
P_TDMA[GRIDAMOUNT-2]=a_E[GRIDAMOUNT-2]/(a_P[GRIDAMOUNT-2]-a_W[GRIDAMOUNT-2]*P_TDMA[GRIDAMOUNT-3]);

Q_TDMA[GRIDAMOUNT-2]=(PHI_L*a_E[GRIDAMOUNT-2]+a_W[GRIDAMOUNT-2]*Q_TDMA[GRIDAMOUNT-3])/(a_P[GRIDAMOUNT-2]-a_W[GRIDAMOUNT-2]*P_TDMA[GRIDAMOUNT-3]);

PHI[GRIDAMOUNT-2]=Q_TDMA[GRIDAMOUNT-2];

for (i=GRIDAMOUNT-3;i>=0;i--)

{
PHI[i]=P_TDMA[i]*PHI[i+1]+Q_TDMA[i];

}

/* 注意：以下为精确解的计算，仅当需要精确解时方才去掉注释语句  */
if (PRECISE){
for (i=0;i<GRIDAMOUNT-1;i++)

{
PHI[i]=(exp(P_Delta*(i+1))-1)/(exp(P_Delta*GRIDAMOUNT)-1)*(PHI_L-PHI_0)+PHI_0;

}
}

//打印输出，并与精确解比较
//  printf("\n PHI\tPHI_PRICISE\tP_TDMA\tQ_TDMA"); //标题行
if((dataoutput=fopen(outputfilename,"wt"))==NULL)

{
printf("\nerror on open output file!");
exit(1);

}
fprintf(dataoutput,"\n %d\t%f",1, (PHI_0-PHI_0)/(PHI_L-PHI_0));//输出相对值 (PHI-PHI_0)/(PHI_L-PHI_0)

for (i=0;i<GRIDAMOUNT-1;i++)

{
fprintf(dataoutput,"\n %d\t%f", i+2, (PHI[i]-PHI_0)/(PHI_L-PHI_0));

}
fprintf(dataoutput,"\n %d\t%f\n",GRIDAMOUNT+1,(PHI_L-PHI_0)/(PHI_L-PHI_0));

fclose(dataoutput);
/*测试strcpy
printf(outputfilename);
printf("\n");
printf(scheme);
printf("\n");
strcpy(outputfilename,scheme);
printf(outputfilename);
printf("\n");
printf(scheme);
printf("\n");
*/

return 0;
}

double max (double number1,double number2)

{
double maxvalue=0;
maxvalue=number2;
if (number1>=number2)

{
maxvalue=number1;
}
return maxvalue;

}

void parseParameter(int *argc, char **argv, double *PHI_0, double *PHI_L, double *P_Delta, char *scheme, char *outputfilename)

{
int i=0;
int pointerofequ=0;

if (*argc<=1)
{
fprintf(stderr,"\n Usage:\n%s –scheme SCHEME\n\t –pdelta P_Delta\n\t –left VALUE-OF-LEFT\n\t –right VALUE-OF_RIGHT\n\t –output OUTPUTFILE\n",argv[0]);

exit(0);
}
for (i=1;i<*argc;i++)

{
if ((argv[i][0]=='-')&&(argv[i][1]=='-'))

{
argv[i]+=2;
pointerofequ=strcspn(argv[i],"=");

if(!strncmp(argv[i],"pdelta",pointerofequ))

{
argv[i]+=pointerofequ+1;
*P_Delta=atol(argv[i]);

}
else if(!strncmp(argv[i],"left",pointerofequ))

{
argv[i]+=pointerofequ+1;
*PHI_0=atol(argv[i]);

}
else if(!strncmp(argv[i],"right",pointerofequ))

{
argv[i]+=pointerofequ+1;
*PHI_L=atol(argv[i]);

}
else if(!strncmp(argv[i],"scheme",pointerofequ))

{
argv[i]+=pointerofequ+1;
strcpy(scheme,argv[i]);

}
else if(!strncmp(argv[i],"output",pointerofequ))

{
argv[i]+=pointerofequ+1;
strcpy(outputfilename,argv[i]);

}
else
{
fprintf(stderr,"\n Unrecognized parameters: –%s\n",argv[i]);

fprintf(stderr,"\n Usage:\n%s –scheme SCHEME\n\t –pdelta P_Delta\n\t –left VALUE-OF-LEFT\n\t –right VALUE-OF_RIGHT\n\t –output OUTPUTFILE\n",argv[0]);

exit(0);
}
}
}
return;

}

• 本文已收录于以下专栏：

举报原因： 您举报文章：一维对流扩散问题TDMA算法 色情 政治 抄袭 广告 招聘 骂人 其他 (最多只允许输入30个字)