为啥数值解与 C 语言中的解析解相同？

Posted 2023-02-18

【中文标题】为啥数值解与 C 语言中的解析解相同？

【英文标题】：Why is the numerical solution coming same as analytical solution in C language?为什么数值解与 C 语言中的解析解相同？

【发布时间】：2019-02-16 10:52:33

【问题描述】：

我编写了一维 cfd 问题，但我的数值解与解析解相同（最多 6 位小数）。

我使用 TDMA 方法进行数值解，而对于解析解，我直接将 x 值代入函数 T(x)。

解析解 T(x) 为T(x) = -(x^2)/2 +11/21(x);

E. G。那么4个网格点；

  x0  = 0.000000,    x1  = 0.333333 ,   x2  = 0.666666 ,   x3  = 0.999999 .
T(x0) = 0.000000 , T(x1) = 0.119048 , T(x2) = 0.126984 , T(x3) = 0.023810.

对于我使用 TDMA 技术的数值解，请参考下面的代码。

输入 n = 4 作为结果。

#include<stdio.h>

void temp_matrix(int n, double *a, double *b, double *c, double *d, double *T);

int main() 
  int Bi = 20.0;
  int n;
  printf("%s ", "Enter the Number of total Grid Points");
  scanf("%d", &n);
  float t = (n - 1);
  double dx = 1.0 / t;
  int i;
  printf("\n");

  double q; // analytical solution below
  double z[n];
  for (i = 0; i <= n - 1; i++) 
    q = (dx) * i;
    z[i] = -(q * q) / 2 + q * (11.0 / 21);
    printf("\nT analytical %lf ", z[i]);
  

  double b[n - 1];
  b[n - 2] = -2.0 * Bi * dx - 2.0;
  for (i = 0; i <= n - 3; i++) 
    b[i] = -2.0;
  

  double a[n - 1];
  a[n - 2] = 2.0;
  a[0] = 0;
  for (i = 1; i < n - 2; i++) 
    a[i] = 1.0;
  

  double c[n - 1];
  for (i = 0; i <= n - 2; i++) 
    c[i] = 1.0;
  

  double d[n - 1];
  for (i = 0; i <= n - 2; i++) 
    d[i] = -(dx * dx);
  

  double T[n];
  temp_matrix(n, a, b, c, d, T);

  return 0;


void temp_matrix(int n, double *a, double *b, double *c, double *d, double *T) 
  int i;
  double beta[n - 1];
  double gama[n - 1];
  beta[0] = b[0];
  gama[0] = d[0] / beta[0];
  for (i = 1; i <= n - 2; i++) 
    beta[i] = b[i] - a[i] * (c[i - 1] / beta[i - 1]);
    gama[i] = (d[i] - a[i] * gama[i - 1]) / beta[i];
  
  int loop;
  for (loop = 0; loop < n - 1; loop++)
    for (loop = 0; loop < n - 1; loop++)

      T[0] = 0;
  T[n - 1] = gama[n - 2];

  for (i = n - 2; i >= 1; i--) 
    T[i] = gama[i - 1] - (c[i - 1] * (T[i + 1])) / beta[i - 1];
  
  printf("\n");
  for (i = 0; i < n; i++) 
    printf("\nT numerical %lf", T[i]);

  

【问题讨论】：

'但我的数值解与解析解相同'.....所以.....这是个问题吗？

旁白：for(i=0;i<=n-1;i++) 风格不佳有两个原因。 for(i=0;i<n;i++) 是惯用的，因为它更具可读性并且对大脑的工作更少，如果n 是无符号的并且0 结束条件将失败。

请正确格式化和缩进您的代码。现在很难读。

为什么数值解在任何给定精度下都与解析解完全匹配是一个问题？当然，人们希望最大限度地提高数值算法重现精确解析解的精度，是吗？

预期结果是什么？

【参考方案1】：

为什么数值解和C语言解析解一样？

它们相差大约 3 位。

以足够的精度打印以查看差异。

使用下面的代码，我们看到 x620 与 T[3] 的 x619 的有效位的最后一个十六进制数字存在差异。这只是 10¹⁵ 差异中的一部分。

#include<float.h>
printf("T analytical %.*e\t%a\n", DBL_DECIMAL_DIG - 1, z[i], z[i]);
printf("T numerical  %.*e\t%a\n", DBL_DECIMAL_DIG - 1, T[i], T[i]);

C 允许在FLT_EVAL_METHOD == 2 时在long double 数学上执行double 数学，然后得到相同的分析/数值结果。由于这一点以及其他细微的 FP 细微差别，您的结果可能与我的不同。

printf("FLT_EVAL_METHOD %d\n", FLT_EVAL_METHOD);

输出

T analytical 0.0000000000000000e+00 0x0p+0
T analytical 1.1904761904761907e-01 0x1.e79e79e79e7ap-4
T analytical 1.2698412698412700e-01 0x1.0410410410411p-3
T analytical 2.3809523809523836e-02 0x1.861861861862p-6

T numerical  0.0000000000000000e+00 0x0p+0
T numerical  1.1904761904761904e-01 0x1.e79e79e79e79ep-4
T numerical  1.2698412698412698e-01 0x1.041041041041p-3
T numerical  2.3809523809523812e-02 0x1.8618618618619p-6

FLT_EVAL_METHOD 0

【讨论】：

@Anshuman Sinha 随着打印精度的提高，您的输出比较如何？

你上面的代码 printf("T 解析 %.*e\t%a\n", DBL_DECIMAL_DIG - 1, z[i], z[i]); 。在我的程序中不起作用。出现以下错误： main.c: In function 'main': main.c:26:43: error: 'DBL_DECIMAL_DIG' undeclared (first use in this function) printf("T analyze %.*e\t%a \n", DBL_DECIMAL_DIG - 1, z[i], z[i]);

@AnshumanSinha 您是否包含#include<float.h>？您使用的是什么 C 编译器/版本？是否符合 C99 标准？ （DBL_DECIMAL_DIG 是在 C99 中引入的）如果仍然有问题，请将DBL_DECIMAL_DIG 替换为17 或(DBL_DIG + 3)

是的，我确实包含了 #include 。我正在使用来自 onlinegdb.com 网站的在线编译器。

运行17.如何解释结果？ T解析1.19047619047619041e-01 0x1.e79e79e79e79ep-4，1.19047619047619041e-01和0x1.e79e79e79e79ep-4一起是什么意思？