Mam prosty program do napisania, kalkuluje szereg fouriera, pakuje współczynniki w array i drukuje dane do stdout, które potem mają być wyświetlone w xmgrace. Współczynniki obliczane są metodą simpsona.
to co na razie napisałem działa dobrze, poza tym, że funkcja obliczona z fouriera ma jakby mniejszą amplitudę.
Nie mogę znaleźć błędu, byłbym bardzo wdzięczny za wszelką pomoc, wskazówki.
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#define PI 3.14159265358979323846
#define PI2 6.28318530717958647692
#define PI_2 1.57079632679489661923
#define INVPI 0.31830988618379067153
double f1 (double x, int k)
{
return x*sin(x*x);
}
double f2 (double x, int k)
{
return f1(x, k)*cos(k*x);
}
double f3 (double x, int k)
{
return f1(x, k)*sin(k*x);
}
/* Metoda Simpsona*/
double integrate( double (*func)(double, int), int k, int n)
{
int i;
static double itg;
double h;
h=(PI2/n);
for( i=1, itg=0 ; i<=n ; i++)
{
itg+=( (h/6)*(func(-PI+i*h,k) + func(-PI+(i-1)*h,k) + 4*func(0.5 * (-PI + i*h + -PI + (i-1)*h),k)) );
}
return itg;
}
/* Obliczanie wspolczynnikow*/
double fourier ( double (*func1)(double, int), double (*func2)(double, int), double (*func3)(double, int), int n, int nmax, double *a, double *b)
{
int k;
b[0]=0;
a[0]=INVPI * integrate ( func1 , 1 , n );
for (k=1 ; k<=nmax ; k++)
{
a[k]=INVPI * integrate ( func2 , k , n );
b[k]=INVPI * integrate ( func3 , k , n );
}
return *a, *b;
}
void fourier_print ( double (*func)(double, int), int nmax, double *a, double *b, int ndiv )
{
int i, n;
static double frr;
for(i=0 ; i<=ndiv ; i++)
{
fprintf( stdout, "%2.6lf %2.6lf\n", (-PI + i*(PI2/ndiv)) , func(-PI + i*(PI2/ndiv), 0) );
}
fprintf(stdout, "&\n");
for(i=0 ; i<=ndiv ; i++)
{
for(n=1, frr=(a[0]/2) ; n<=nmax ; n++)
{
frr+=( a[n]*cos(n*(-PI + i*(PI2/ndiv))) + b[n]*sin(n*(-PI + i*(PI2/ndiv))) );
}
fprintf( stdout, "%2.6lf %2.6lf\n", (-PI + i*(PI2/ndiv)) , frr );
}
}
int main ( int argc, char *argv[] )
{
int nmax, ndiv, n;
double *a;
double *b;
fprintf(stderr, "Enter n: "); scanf("%i", &n);
fprintf(stderr, "Enter nmax: "); scanf("%i", &nmax);
fprintf(stderr, "Enter ndiv: "); scanf("%i", &ndiv);
a = malloc ( (nmax+1) * sizeof(double) );
b = malloc ( (nmax+1) * sizeof(double) );
fourier ( f1, f2, f3, n, nmax, a, b);
fourier_print ( f1 , nmax, a, b, ndiv );
return EXIT_SUCCESS;
}