一维对流扩散问题TDMA算法

转载 2012年03月21日 22:05:10
/*
一维对流扩散问题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
参数顺序无关紧要。

注意现在的代码输出相对值(phi-phi0)/(phi_L-phi0),想输出绝对值phi的自己改代码吧:)
另存为upwind.c编译即可
//*/
#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;
	
}

算法/贪心算法/BinPacking一维装箱问题

问题描述We have a number of bins each with a capacity of 1, and we have a set of objects all with differ...
  • Carinya
  • Carinya
  • 2017年04月25日 20:00
  • 1458

装箱问题近似算法概述

问题描述:一维经典装箱问题可描述如下:S=(S1,S2,..Sn),其中0< Si ≤ 1, 称之为第i个物体的体积(或重量),1≤i≤n,现有n个容积(或载重量)为1 的箱子,要求如何设法将S1,S...
  • x_i_y_u_e
  • x_i_y_u_e
  • 2015年07月05日 19:52
  • 2710

分治算法 最接近点对(一维)

最接近点对,在实际应用中有很强的实用性!下面我就介绍一下,怎么求一维的最接近点!二维的会再下一次博客更新! 方法:把所有点从大到小先进行排序,找出最大点MAX和最小点min。然后求mid=(...
  • hengliwuyou
  • hengliwuyou
  • 2015年06月19日 16:15
  • 1862

一维数组解决01背包问题

01背包是在M件物品取出若干件放在空间为W的背包里,每件物品的体积为C1,C2,…,Cn,与之相对应的价值为W1,W2,…,Wn.求解将那些物品装入背包可使总价值最大。 设置一个数组,V[1...
  • amanicspater
  • amanicspater
  • 2015年11月19日 22:25
  • 816

三种多址技术FDMA CDMA TDMA 的理解

TD-SCDMA、WCDMA、CDMA2000三种方式的共性是多址方式用都用到CDMA,三者最主要的区别是双工方式不同,TD-SCDMA是时分双工TDD,WCDMA和CDMA2000是频分双工FDD。...
  • haijin0829
  • haijin0829
  • 2016年03月24日 12:27
  • 3362

一点一滴完全突破KAZE特征检测算法,从各向异性扩散滤波开始(1)

ECCV2012中出现了一种比SIFT更稳定的特征检测算法KAZE。尽管,这个算法是几个法国人提出的,但是算法却有一个日文的名字。KAZE是日语‘风’的谐音,最近宣布退休的宫崎骏所拍摄的影片“起风了”...
  • baimafujinji
  • baimafujinji
  • 2013年09月11日 17:31
  • 8602

编程珠玑之第二章questionB: n元一维向量旋转问题

问题描述: B.将一个n元一维向量向左旋转i个位置。例如,当n=8且i=3时,向量abcdefgh旋转为defghabc. 简单的代码使用一个n元的中间向量在n步内完成该工作。你能否仅使用数十个额外字...
  • JohnnyHu90
  • JohnnyHu90
  • 2015年01月07日 14:56
  • 975

MAC TDMA网络的实现

MACTDMA系统的设计与实现         为实现传输的实时性,设计并实现了一个小型的TDMA系统。...
  • shenlei190810
  • shenlei190810
  • 2014年12月05日 12:24
  • 1806

各向异性扩散PM模型原理与C++实现

本文介绍了各向异性扩散PM模型,并给出了C++代码实现。 一、PM模型原理 其中,                                     ...
  • cyh706510441
  • cyh706510441
  • 2015年04月24日 16:22
  • 3060

NS2.29中Tdma的实现分析

1.  NS2中的Tdma时帧结构如下图: 其中前导的数据结构为:static int *tdma_preamble_; tdma_preamble_ = new int[max_slot_n...
  • u010693479
  • u010693479
  • 2013年08月19日 10:48
  • 1259
内容举报
返回顶部
收藏助手
不良信息举报
您举报文章:一维对流扩散问题TDMA算法
举报原因:
原因补充:

(最多只允许输入30个字)