我正在尝试编写从-inf到inf的Romberg积分的c代码。这里写的是我给予数字并运行函数的地方。当大于20的数字是15或有时也是15时，代码不会运行或花费太长时间。
我很确定这是因为功能，但我真的不知道哪个部分应该改变。
先谢谢你了
#include <stdio.h>
#include <stdlib.h>
#include <limits.h>
#include <math.h>
#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); //HERE
fprintf(fp,"%lf\t%25.16e\n",z,r*alpha/a[i]);
}fclose(fp);
}
// free(R);
return 0;
}我将变量改为 var而不是var，并且只使用两个一维数组而不是R[n][n]二维数组。
++将Rpre和Rcur改为R，inf改为2^27，for-循环的int x = 1 << (i - 1); for (int j = 1; j <= x; j++) { pow（x，2）改为xx
2条答案
按热度按时间szqfcxe21#
我很确定这是因为功能，但我真的不知道哪个部分应该改变。
不，它是您正在使用的算法的特征。相对于
romberg()的最后一个参数的值，您的计算扩展得非常快（O（2n）），这就是您所询问的。对于任何给定的函数，在程序运行的时间超过您愿意等待的时间之前，n 不必非常大，这并不奇怪。你可以通过代码和编译的改进来稍微增加这个限制，但是即使使它的速度提高一倍，你也只能在任何固定时间内计算的 n 中增加1。
您还应该考虑一些计算方面的考虑，尤其是浮点表示。FP表示法并不精确，而且您选择的积分界限几乎肯定是两者
同时，使用如此大的边界意味着你需要去相当大的
n来获得积分的良好估计，因为它在整个范围中的比例很小，对该估计有意义。可能的算法改进包括：
double有53位尾数（大多数实现都是这样），这会让你刚好超过你开始失去精度的边缘。还有其他选择。ccgok5k52#
一些调整（产生1.7倍加速）...
1.执行
(x * x)比pow(x,2)快1.对整数执行
pow(2,j)是很慢的--它可以用移位来代替。1.它应该只做 * 一次 *（并且 * 不是 * 在每次循环迭代中）。
1.在每次循环迭代中执行
n_k = pow(4,j)很慢1.可以替换为
n_k *= 4.0这是正确的代码。它标注了bug和修复。我还在
printf调用中添加了基准计时：在上面的代码中，我使用了
cpp条件来表示old与新代码：还有...
注意：这可以通过
unifdef -k -DORIG=0运行文件来清理。循环迭代的原始版本和固定版本的运行时间：
加速比为1.699x
下面是原始程序 * 和 * 修复程序（用
-DDEMO=10编译）的输出。