yxl1.cpp

用差分方法求解沃尔泰拉方程

/*Solveing differential equation
  dy/dt=Lamda*Integrate[(t-s)^(Alpha)*y(s),{s,0,t}]*/

#include <stdlib.h>
#include <stdio.h>
#include <values.h>
#include <dos.h>
#include <time.h>
#include <math.h>
#include <fcntl.h>
#include <stat.h>
#include <string.h>
#include<io.h>


#define LAMBDA  	-1
#define ALPHA    	1
#define MAXN     	5000
#define DEFAUTN  	100
#define T0			0
#define T1			1
#define Y0			1


double   y[MAXN];
double increase=0,h;
  int maxn=DEFAUTN;
  double lambda=double(LAMBDA),
	 alpha=double(ALPHA),
	 t0=double(T0),
	 t1=double(T1),
	 y0=double(Y0);
  int n,k;
  char filename[80]={"\\yxl1"};
  FILE * fp;
  struct time tm;
  struct date dt;



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


  if(argc==2){
    strcat(filename,argv[1]);
    maxn=atoi(argv[1]);
  }
  else if(argc>2){
    printf("command line arguments error!\n");
    return 1;
  }
  if(maxn>=MAXN){
    printf("Input n is too large!\n");
    return 1;
  }
  strcat(filename,".txt");

  if((fp = fopen(filename,"w+t"))==NULL){
    printf("open file error!");
    return 1;
  }

  getdate(&dt);
  if(fprintf(fp,"Current date:\t %d.%d.%d\n",
		    dt.da_year,dt.da_mon,dt.da_day)==EOF){
    printf("write file error");
    return 2;
  }
  gettime(&tm);
  if(fprintf(fp,"Current time:\t%2u:%02u:%02u\n",
		 tm.ti_hour,tm.ti_min,tm.ti_sec)==EOF){
    printf("write file error!");
    return 2;
  };
  if(fprintf(fp,"n=%d alp=%f lam=%f T0=%f T1=%f Y0=%f\n",
	      maxn,alpha,lambda,t0,t1,y0)==EOF){
    printf("write file error!");
    return 2;
  };
  if(fprintf(fp, "%8s %15s %15s\n",
	    "step","value","error")==EOF){
    printf("write file error!");
    return 2;
  };


  y[0]=y0;
  h=double((t1-t0))/double(maxn);
  for(n=1;n<=maxn;n++){
     increase=0;
    for(k=0;k<n;k++)
      increase+=double(n-k)*h*y[k];
    y[n]=y[n-1]+lambda*h*h*increase;
    if(fprintf(fp, "%8d %15f\n",n,y[n])==EOF){
	   printf("write file error!");
      return 2;
    }
  }
  fclose(fp);
  return 0;
}