c - C 语言中的一维自适应网格细化

Question

我是一名物理学博士生，目前正在进入数值相对论领域，我必须在二维中实现自适应网格细化代码。与程序的其他部分一样，我通常更喜欢做一些更简单的事情来了解正在发生的事情，然后再跳到更复杂的情况。然而，我似乎仍然在做一些根本性错误的事情。

我的代码执行(或至少应该执行)以下过程:我以 N 个大小为 h 的间隔离散化 x 轴。每次计算一个点时，程序都会停止并通过将大小为 h 的间隔更改为另一个具有两个步长 h/2 的间隔来再次计算该点。该程序检查结果是否低于某个用户指定的容差，如果没有，则该过程以步长 h/4 再次开始，依此类推。下面的草图说明了该过程

细化函数起作用后，我绝对没有兴趣将函数的值保留在细化网格上。我想要的只是以最高精度计算粗网格上的函数(在图像中我想要保留和更改的是粗[基本]网格的黑点的值)。

不幸的是，在通过细化算法后，我发现解决方案没有任何改进。我不期望函数的绘图是完美的，但我期望每个点都非常接近解析解。这是我的细化函数(该函数被递归调用，直到达到用户指定的最大细化级别):

void refine( int l, long double dx, long double x_min, long double x_max, long double f_min, long double *f_max ){
  // l = level of refinement, dx = step size, x_min is current x position, x_max = point we want to calculate, f_min = function evaluated at x_min, f_max = function evaluated at x_max
  int i;
  long double *f_aux, f_point;

  f_aux = (long double *) malloc ( (2*l + 1) * sizeof (long double) );

  dx = 0.5 * dx;

  f_aux[0] = f_min;
  for( i=1; i<2*l+1; i++ ){
    f_aux[i] = ( 1.0 - 2.0 * dx * ( x_min + (i-1)*dx - X0 ) / DELTA ) * f_aux[i-1];
  }

  if( l < lMAX ){
    if( fabs( f_aux[2*l] - *f_max ) > TOL ){
      f_point = f_aux[2*l];
      free( f_aux );
      l++;
      refine( l, dx, x_min+dx, x_max, f_min, &f_point );
    }
    else{
      *f_max = f_aux[2*l];
      free( f_aux );
    }
  }
  else{
    *f_max = f_aux[2*l];
    free( f_aux );
  }

  return;
}

有人能解释一下这个问题吗？我感觉完全被困住了。

提前致谢!

Answer 1

看起来您的迭代在最后一点附近有所增强，但在起点附近仍然很粗糙:

*---------------* refine 0

*-------*-------* refine 1

*-------*---*---* refine 2

*-------*---*-*-* refine 3

由于您的方程看起来是双曲线的(取决于之前的所有迭代)，因此解误差本质上将是累积的。我认为你应该从一开始就使用精细网格进行迭代以获得正确的解决方案 - 不仅仅是在附近。例如。对于具有此类实现的简单方程 df/dx=f:

float update(float f_prev, float dx, float /* unused */ target_dx) {
    return f_prev*dx + f_prev;
}

int main(void)
{
    int n = 8;
    float dx = 2.f/(float)(n-1);
    float f[n];
    f[0] = 1.f;

    for (int i = 1; i < n; i++) {
        f[i] = update(f[i-1], dx, dx / 16.f);
    }

    return 0;
}

我会使用简单的递归公式:

float update_recursive(float f_prev, float dx, float target_dx) {
    if (dx > target_dx) {
        return update_recursive(
            update_recursive(f_prev, dx / 2.f, target_dx), dx / 2.f, target_dx
        );
    }

    return dx * f_prev + f_prev;
}

该方法提高了解决方案的质量。终止条件可能比 target_dx 更适合解决方案。当然，需要保证递归是有界的。