傅里叶级数系数gnuplot图

Question

这个程序的目的是找出傅里叶级数的系数，然后根据计算值绘制函数。

分片是

h(t)=1+(t/pi)表示t<0 h(t)=1-(t/pi)表示t>0

#include <stdio.h>
#include <math.h>
#include "comphys.c"
#include "comphys.h"
#define pi 3.141592653589793

double trap1 (double l, double u, int m, int N);
double trap(double l,double u,int m,int N); // Use trapezoidal method
double inte(double x,int m);
double integrand(double x, int m); // This is the integrand of the bms
double w=1.0; //omega (here - w-2pi/Period where Period=2pi)
int main()
{
  FILE *out,*rsp;
  int m; // Index of the b coefficients
  int M=5; // Number of coeeficients to compute
  int T=201;  // Number of time steps 
  int t; // This is a counter to keep track of the h(t) function to plot
  double *b,*h,j,dj,*a,a0;
  double ll,ul;
  a=dvector(1,M);
  b=dvector(1,M);  // coefficients
  h=dvector(1,T);  // function built from Fourier terms
  ll=-pi;  // lower time limit
  ul=pi;   // upper time limit

  // Prepare to write for plotting
  printf("\nBrute force complex Fourier transform algorithm");
  if ((out=fopen("fs.out","w"))==NULL) {
    printf("\nCannot open file for output\n");
    getchar();
    return(0);          
  }

  // Loop through the wave numbers    
  for (m=1;m<=M;m++)
    {
      b[m]=(1.0/pi)*trap(ll,ul,m,T);
      a[m]=(1.0/pi)*trap1(ll,ul,m,T);
      printf("\nb[%d]=%f   a[%d]=%f",m,b[m],m,a[m]);

      // Build the h(t) function as we go
      t=0;
      dj=(ul-ll)/(double)T;
      for (j=ll;j<=ul;j+=dj)
    {
      t++;
      h[t]+=(b[m]*sin((double)m*w*j))+(a[m]*cos((double)m*w*j));
    }
    }

  // Write out for plotting
  t=0;
  for (j=ll;j<=ul;j+=dj)
    {
      t++;
      fprintf(out,"%f %f\n",j,h[t]);
    }
  fclose(out);
  printf("\nProgram complete without known error.\n");

  if((rsp = fopen("gnuplot.rsp","w")) == NULL)
    {
      printf("\nCannot open gnuplot.rsp for writing\n");
      exit(1);
    }
  fprintf(rsp,"plot '%s' using 1:2 with lines\n",out);
  fprintf(rsp,"pause mouse\n");
  fprintf(rsp,"replot\n");
  fclose(rsp);
  if(system("gnuplot gnuplot.rsp") == -1)
    {
      printf("\nCommand could not be executed\n");
      exit(1);
    }
  return;
}

double trap (double l, double u, int m, int N)

{
  double x,integral,dx;
  int i;

  if (u==l) return (0.0);
  dx=(u-l)/(double)(N-1);

  integral=0.0;
  integral+=0.5*(integrand(u,m)-integrand(l,m));
  x=l+dx;
  for (i=1;i<N;i++)
    {
      integral += integrand(x,m);
      x+=dx;
    }
  integral*=dx;
  printf("\n%f\n",integral);
  return(integral);
}

double trap1 (double l, double u, int m, int N)
{
  double x,integral,dx;
  int i;

  if (u==l) return (0.0);
  dx=(u-l)/(double)(N-1);

  integral=0.0;
  integral+=0.5*(inte(u,m)-inte(l,m));
  x=l+dx;
  for (i=1;i<N;i++)
    {
      integral += inte(x,m);
      x+=dx;
    }
  integral*=dx;
  return(integral);
}


double integrand(double x, int m)
{
  if (x<0) return (sin((double)m*w*x)+((x/pi)*sin((double)m*w*x)));
  return(sin((double)m*w*x)-((x/pi)*sin((double)m*w*x)));
}

double inte(double x,int m)
{
  if (x<0) return (cos((double)m*w*x)+((x/pi)*cos((double)m*w*x)));
  return(cos((double)m*w*x)-((x/pi)*cos((double)m*w*x)));
}

我的问题是，它会很好地编译并列出系数(不知道它们是否正确或其他什么)，但在试图绘制它时，它给了我最后的答案。

“无法读取数据文件”、"gnuplot.rsp"，行:1 util.c:没有此类文件或目录

我不知道怎么解决这个问题。我用gnuplot绘制了很多程序，所以它不应该是缺少的文件或任何东西。

Answer 1

这里：

fprintf(rsp,"plot '%s' using 1:2 with lines\n",out);  

您应该直接使用字符串"fs.out“而不是"out”，这是一个文件指针(并且已经关闭)。