计算科学 - 如何通过与 GNU Scientific Library (GSL) 集成来获得乌龟坐标？ - 吾爱随笔录

如何通过与 GNU Scientific Library (GSL) 集成来获得乌龟坐标？

计算科学计算物理学坐标变换 gsl

2021-12-21 21:50:05

我试图获得整合表达式的乌龟坐标（Eddington-Finkelstein Coordinates）的值： $\frac{dr^*}{dr} = (1 - rS/r)^{-1}$ 使用 GNU 科学图书馆 (GSL)。我知道解析表达式是 $r* = r + rs·ln(r/rS - 1)$ 所以我可以比较两个结果：分析的和综合的。

不幸的是，集成解决方案似乎轮换了

但即使我旋转图表，结果也会发生变化……

如果我添加 2 个常量，则图表是相同的

我的问题是：

为什么计算是旋转的？
为什么如果我旋转图表它仍然移动？（值 -x + 10.605, -yxCalculated[i] + rS）两个图重叠）

我需要知道如何进行这种计算，因为在黑洞反弹或虫洞的情况下，我们通常没有解析表达式。

这是我的代码：

#include <stdio.h>
#include <assert.h>
#include <math.h>

#include <gsl/gsl_errno.h>
#include <gsl/gsl_odeiv2.h>

#define SYS_DIM 1
#define NODES 10001
#define SCH_MIN_INF -10000
#define SCH_PLU_INF 10000

struct yxParams
{
  double rS;  /* Schwarzschild radius */
};
double rS = 2;
/*
  I refer to r* as y and r as x:
    r = x
    r* = y
  So yx = r*(r) and xy = r(r*)
*/
double yxAnalytical[NODES];
double yxCalculated[NODES];

/**
  @brief  ODE IV 1st grade to solve y[x]'

  @param  t     Independent variable
  @param  y[]   Left side of the 1st grade system of equations
  @param  sys[] Right side of the 1st grade system of equations
*/
int
funcY(double t,
      const double y[],
      double sys[],
      void* params)
{
  struct yxParams *par = (struct yxParams*)params;
  double rS = (par->rS);

  (void)(t);
  double x0 = y[0];

  assert(x0 != 0);
  sys[0] = 1.0/(1.0 - rS/x0);

  return GSL_SUCCESS;
}

/**
  @brief  r*(r) integration... so the inverse r(r*) tortoise coordinates.
          We do this because it's easier to get an analytic expression from
          r*(r) that from r(r*)

  @param  IN, a: interval lower limit
  @param  IN, b: interval upper limit
  @param  IN, n: number of subintervals
  @param  IN, ic: Initial condition
  @param  OUT, yx[]: Array with the integration values
*/
int
yx_integration(double a, double b, int n, double iC)
{
  int status = GSL_SUCCESS;

  double h = (b - a)/n;

  double x0 = a;
  double x1 = x0 + h;
  double epsAbs = 0;
  double epsRel = 1e-6;

  struct yxParams params = {rS};

  gsl_odeiv2_system sys = {funcY, NULL, SYS_DIM, &params};
  const gsl_odeiv2_step_type* T = gsl_odeiv2_step_rkf45;
  gsl_odeiv2_driver* d = gsl_odeiv2_driver_alloc_y_new (&sys, T, h, epsAbs, epsRel);

  double y[1] = {iC};
  int i;
  for (i = 0; i < NODES; i++)
  {
    status = gsl_odeiv2_driver_apply(d, &x0, x1, y);
    x0 = x1;
    x1 = x0 + h;

    if (status != GSL_SUCCESS)
    {
      printf ("error, return value = %d\n", status);
      break;
    }
    else
    {
      double tmp = y[0];

      if(_isnanf(tmp) || isinf(tmp))
      {
        tmp = SCH_MIN_INF;
      }

      yxCalculated[i] = tmp;
    }
  }
  gsl_odeiv2_driver_free(d);

  return status;
}

/**
  @brief  Tortoise analytical expression: r*(r)

  @param  IN a: interval lower limit
  @param  IN h: step
*/
void
yx_analytical(double a, double h)
{
  double x = a;
  int i;
  for(i = 0; i < NODES; i++)
  {
    double tortoise = x + rS*log(x/rS - 1);

    if(_isnanf(tortoise))
    {
      tortoise = SCH_MIN_INF;
    }
    yxAnalytical[i] = tortoise;
    x += h;
  }
}

/**
  @brief  Save the data so we can plot it

  @param  IN x
  @param  IN h
*/
void
txt_data(double x, double h)
{
  FILE *fp;
  fp = fopen("data.txt", "w");

  fprintf(fp, "#%32s %32s %32s %32s\n", "yxAnalytical", "x", "yxCalculated", "x");
  int i;
  for (i = 0; i < NODES; i++)
  {
    /*
    If we plot 'yx vs x', we're plotting the inverse function... so we're plotting
    r(r*) the tortoise coordinates xy[].

    Q1) yxCalculated is rotated with respect yxAnalytical. Why?
    Q2) if I plot (if I'm not wrong this implies rotating 180º the graph)
        fprintf(fp, "%.32f %.32f %.32f %.32f\n", yxAnalytical[i], x, -x, -yxCalculated[i]);
        both graphs are practically the same... but not equal. I need to shift the graph
        fprintf(fp, "%.32f %.32f %.32f %.32f\n", yxAnalytical[i], x, -x + 10.605, -yxCalculated[i] + rS);
        so calculations seems to be in the right direction... but certainly something is really wrong
    */
    fprintf(fp, "%.32f %.32f %.32f %.32f\n", yxAnalytical[i], x, yxCalculated[i], x);
    /*
      Trick to overlap the graphs...
    *
    fprintf(fp, "%.32f %.32f %.32f %.32f\n", yxAnalytical[i], x, -x + 10.605, -yxCalculated[i] + rS);
    */
    x += h;
  }
  fclose(fp);
}

/**-----------------------------------------------------------------------------
------------------------------------------------------------------------------*/
int main()
{
  /*
  To our purposes these values are enough far away from the BH.
  These values are for the numerical integration.
  */
  double a = -150;
  double b = 150;
  int n = NODES - 1;
  double h = (b - a)/n;

  /* Do calculations... */
  yx_analytical(a, h);
  yx_integration(a, b, n, a);

  /* Save the data in a text file */
  txt_data(a, h);

  return 0;
}

我会自己回答

这是我的最终代码（向后集成）。我添加了 jacobian 并用我的 funcY 解决了一个小问题（误解）。

#include <stdio.h>
#include <assert.h>
#include <math.h>

#include <gsl/gsl_errno.h>
#include <gsl/gsl_odeiv2.h>
#include <gsl/gsl_matrix.h>

#define SYS_DIM 1
#define NODES 10001
#define SCH_MIN_INF -10000
#define SCH_PLU_INF 10000

struct yxParams
{
  double rS;  /* Schwarzschild radius */
};
double rS = 2;
/*
  I refer to r* as y and r as x:
    r = x
    r* = y
  So yx = r*(r) and xy = r(r*)
*/
double yxAnalytical[NODES];
double yxCalculated[NODES];

/**
  @brief  yx[] jacobian, so you can use implicit methods

  rS/(r^2*(1 - rS/r)^2)
*/
int
jacoY(double r,
     const double y[],
     double *dfdy,
     double dfdt[],
     void *params)
{
  (void)(y);
  struct yxParams *par = (struct yxParams*)params;
  double rS = (par->rS);

  gsl_matrix_view dfdy_mat = gsl_matrix_view_array(dfdy, 1, 1);
  gsl_matrix *m = &dfdy_mat.matrix;
  /*gsl_matrix_set (m, 0, 0, -rS/(r*r*(rS/r - 1)*(rS/r - 1)));*/
  gsl_matrix_set (m, 0, 0, -rS/(r*r*(1 - rS/r)*(1 - rS/r)));
  dfdt[0] = 0.0;

  return GSL_SUCCESS;
}

/**
  @brief  ODE IV 1st grade to solve y[x]'

  @param  x     Independent variable
  @param  y[]   Left side of the 1st grade system of equations
  @param  sys[] Right side of the 1st grade system of equations
*/
int
funcY(double x,
      const double y[],
      double sys[],
      void* params)
{
  struct yxParams *par = (struct yxParams*)params;
  double rS = (par->rS);

  assert(x != 0);

  sys[0] = 1.0/(1.0 - rS/x);

  return GSL_SUCCESS;
}

/**
  @brief  r*(r) integration... so the inverse r(r*) tortoise coordinates.
          We do this because it's easier to get an analytic expression from
          r*(r) that from r(r*)

  @param  IN, a: interval lower limit
  @param  IN, b: interval upper limit
  @param  IN, n: number of subintervals
  @param  IN, ic: Initial condition
  @param  OUT, yx[]: Array with the integration values
*/
int
yx_integration(double a, double b, int n, double iC)
{
  int status = GSL_SUCCESS;

  double h = (b - a)/n;

  double x0 = a;
  double x1 = x0 + h;
  double epsAbs = 0;
  double epsRel = 1e-6;

  struct yxParams params = {rS};

  gsl_odeiv2_system sys = {funcY, jacoY, SYS_DIM, &params};
  const gsl_odeiv2_step_type* T = gsl_odeiv2_step_rk8pd;  /* rkf45, rkck, rk8pd; rk4imp bsimp msadams msbdf */
  gsl_odeiv2_driver* d = gsl_odeiv2_driver_alloc_y_new (&sys, T, h, epsAbs, epsRel);

  double y[1] = {iC};
  int i;
  for (i = 0; i < NODES; i++)
  {
    status = gsl_odeiv2_driver_apply(d, &x0, x1, y);
    x0 = x1;
    x1 = x0 + h;

    if (status != GSL_SUCCESS)
    {
      printf ("error, return value = %d\n", status);
      printf("x0: %f, y[0]: %f", x0, y[0]);
      break;
    }
    else
      yxCalculated[i] = y[0];
  }
  gsl_odeiv2_driver_free(d);

  return status;
}

/**
  @brief  Tortoise analytical expression: r*(r)

  @param  IN a: interval lower limit
  @param  IN h: step
*/
void
yx_analytical(double a, double h)
{
  double x = a;
  int i;
  for(i = 0; i < NODES; i++)
  {
    double tortoise = x + rS*log(x/rS - 1);

    if(isnan(tortoise))
      tortoise = SCH_MIN_INF;
    yxAnalytical[i] = tortoise;
    x += h;
  }
}

/**
  @brief  Save the data so we can plot it

  @param  IN x
  @param  IN h
*/
void
txt_data(double x, double h)
{
  FILE *fp;
  fp = fopen("data.txt", "w");

  int n = NODES - 1;
  int i;
  for (i = 0; i < NODES; i++)
  {
    fprintf(fp, "%.32f %.32f %.32f %.32f\n", yxAnalytical[i], x, yxCalculated[n - i], x);
    x += h;
  }
  fclose(fp);
}

/**-----------------------------------------------------------------------------
------------------------------------------------------------------------------*/
int main()
{
  /*
  To our purposes these values are enough far away from the BH.
  These values are for the numerical integration.
  */
  double a = -150;
  double b = 150;
  int n = NODES - 1;
  double h = (b - a)/n;

  /* Do calculations... */
  yx_analytical(a, h);
  /*yx_integration(b, a, n, b);*/
  /* Unfortunately we need and exact IC... we cannot use the tortoise approximation
  rS = r enough far away... */
  double ic = b + rS*log(b/rS - 1);
  yx_integration(b, a, n, ic);
  /* Save the data in a text file */
  txt_data(a, h);

  return 0;
}

因此，如果我们对初始条件使用乌龟近似值，其中 rS = r 很远，我们就会得到一条具有完美形状的曲线......但是有一个偏移注意，在地平线上我们不能继续积分（但是好的）。如果我们使用确切的 IC，其中 $r_* = r + rS·ln(r/rS - 1)$ 我们得到了准确的曲线，除了我们不能在地平线上继续整合。

1个回答

理论上你的初始条件到底是什么？为什么在评估否定参数的对数时不会出现运行时错误？正如您在代码中强加的那样 $y(a)=a$ 为了 $a=-150$ ，你当然会得到一个解决方案 $y(x)=x+\log((x-r_s)/(a-r_s))$ （见下文），这明显不同于在 $x=a$ .

再说一次，你实际求解的微分方程是

y^{'} (x) = \frac{1}{1 - r_{s} / y (x)} = \frac{y (x)}{y (x) - r_{s}} ⟹ y^{'} (x) - r_{s} \frac{y^{'} (x)}{y (x)} = 1 ⟹ y (x) - r_{s} \ln | y (x) | = x + c ⟹ y (x) e^{- y (x) / r_{S}} = C e^{- x / r_{S}}

$y'(x)=\frac1{1-r_s/y(x)}=\frac{y(x)}{y(x)-r_s} \\ \implies y'(x)-r_s\frac{y'(x)}{y(x)}=1 \\ \implies y(x)-r_s\ln|y(x)|=x+c \\ \implies y(x)e^{-y(x)/r_S}=Ce^{-x/r_S}$ 这可以使用 Lambert-W 函数来解决。（如果

y (a) = a

$y(a)=a$ ，然后

C = a

$C=a$ 跟随。）

为了提供对比，从问题文本中猜测的 ODE 或积分问题读作

\frac{d y (x)}{d x} = \frac{1}{1 - r_{S} / x} = \frac{x}{x - r_{S}} = 1 + \frac{r_{S}}{x - r_{S}} .

$\frac{dy(x)}{dx}=\frac1{1-r_S/x}=\frac{x}{x-r_S}=1+\frac{r_S}{x-r_S}.$ 这很好地集成到

r^{*} = y (x) = x + r_{s} \ln | x - r_{S} | + C

$r^*=y(x)=x+r_s\ln|x-r_S|+C$ 无论初始条件是什么，它都将参考解决方案的公式作为其案例之一。

要实现这一点，您需要注意什么是因变量和自变量。那么 ODE 函数应该被实现为

int
funcY(double x,
      const double y[],
      double sys[],
      void* params)
{
  struct yxParams *par = (struct yxParams*)params;
  double rS = (par->rS);

  assert(x != rS);
  sys[0] = x/(x - rS);

  return GSL_SUCCESS;
}

其它你可能感兴趣的问题

上一篇数值偏微分检查下一篇比较两个不同函数的最小值