我正在尝试为从 -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] 二维数组。