C 集成运行太慢

问题描述 投票:0回答:0

我正在尝试为从 -inf 到 inf 的 Romberg 集成编写 C 代码。 “这里”写的是我给出数字并运行函数的地方。 当数字超过 20 或有时也有 15 时,代码不会运行或运行时间太长。

我很确定这是因为功能,但我真的不知道应该更改哪个部分。

提前致谢

/**************************************
*kompile
*gcc -Wall -pedantic h21.c -o h21 -lm
*aufrufen mit
*./h21
*******************************************/
#include <stdio.h>
#include <stdlib.h>
#include <limits.h>
#include <math.h>
//#define N LONG_MAX
#define eps pow(10,-8)
double f (double a, double x, double z){
    return exp(-pow(x,2)/pow(a,2))/sqrt(pow(x,2)+pow(z,2));
}
                       
long double romberg(double f(double, double, double), double a, double xa, double xe, double z,int N){

    double *Rpre=malloc(N*sizeof(double));
    double *Rcur=malloc(N*sizeof(double));
    //double R1[n],R2[n];
    //double Rpre=&R1[0],*Rcur=&R2[0]; //Rpre fuer previous row, Rcur fuer current row
    double h =xe-xa;
    Rpre[0]=0.5*h*(f(a,xa,z)+f(a,xe,z)); //erste Trapaze
    
    for(int i=1;i<N;i++){
        h*=0.5;
        double sum=0;
        
        for(int j=1;j<=pow(2,i-1);j++){//bis 2^(n-1)
            sum+=f(a,xa+(2*j-1)*h,z);
        }
        Rcur[0]=h*sum+0.5*Rpre[0];//R[i][0]
        
        for(int j=1;j<=i;j++){
            double n_k=pow(4,j);
            Rcur[j]=(n_k*Rcur[j-1]-Rpre[j-1]) /(n_k-1); //R[i][j]
        }
        
       /* if(i>1&&fabs(Rpre[i-1]-Rcur[i])<eps){ //Wenn die rel. Genauigkeit erfuellt ist.
            return Rcur[i];
        }*/
       double *temp=Rpre; //wenn es bis zum Ende der Reihe laeuft, muss Current Wert Previous Wert sein
        Rpre=Rcur;
        Rcur=temp;
        
    }printf("%25.16e\n",Rpre[N-1]);
    return Rpre[N-1];
            

        
}
            
int main(int argc, char **argv) {
    double a[4]={0.02,1,2,4};
    double N=LONG_MAX; //maximale double wert
    double alpha=1/sqrt(M_PI);
    
    //long double ** R= malloc( N * sizeof(long double ) );
//printf("%lf",N);
    
 /* Parameter für den z-Bereich 
   */
  double const zmin = 0.;
  double const zmax = 5.;
  int const znum = 100;

  /* Schleife über z-Werte, für die integriert werden soll
   */
for (int i=0;i<4;i++){
    char filename[32];
    sprintf(filename, "h21_a_%.2lf.txt",a[i]);
    FILE *fp = fopen(filename, "w");
    
 for ( int iz = 0; iz <= znum; iz++ )
  {

    double z = zmin + ( zmax - zmin ) / (double)znum * iz ;
     double r=romberg(f,a[i],-N,N,z,20);
    fprintf(fp,"%lf\t%25.16e\n",z,r*alpha/a[i]);
     

  }fclose(fp);
}
   // free(R);
    return 0;
}
            

我将变量更改为 *var 而不是 var 并且只使用两个一维数组而不是 R[n][n] 二维数组。

c runtime integration
© www.soinside.com 2019 - 2024. All rights reserved.