为什么 C 为这个算法赋予 Python 不同的值？

Question

我用 C 编写了一个简短的程序来执行线性插值，它迭代以将函数的根赋予给定的小数点数。到目前为止，对于函数 f(x):

    long double f(long double x) {
        return (pow(x, 3) + 2 * x - 8);
    }

程序无法收敛到 1dp 值。程序更新变量 a 和 b，f(x) 的根位于它们之间，直到 a 和 b 都以给定的精度舍入到相同的数字。使用 long double 和上述函数，调试器显示前 2 次迭代：

    a = 1.5555555555555556
    a = 1.6444444444444444

虽然应该是：

    a = 1.5555555555555556
    a = 1.653104925053533

此后程序无法更新值。我使用的线性插值方程是给定 here 数学方程的重新排列版本，我使用的代码是我编写的 python 程序的 C 版本。为什么 C 实现得到不同的值，尽管算法相同，我该如何解决？

好的，我仍然掌握了这个窍门，但希望下面我有一个最小的、完整的和可验证的示例：

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

long double a; long double b; long double c; // The values for interpolation
long double fofa; long double fofb; long double fofc; // The values f(a), f(b) and f(c)
const int dp = 1; // The number of decimal places to be accurate to

long double f(long double x) {
    return (pow(x, 3) + 2 * x - 8);
}

int main(void) {
    a = 1; b = 2;
    while(roundf(a * pow(10, dp))/pow(10, dp) != roundf(b * pow(10, dp))/pow(10, dp)) { // While a and b don't round to the same number...

        fofa = f(a); fofb = f(b); // Resolve the functions
        printf("So f(a) = %g, f(b) = %g\n", (double)fofa, (double)fofb); // Print the results

        c = (b * abs(fofa) + a * abs(fofb)) / (abs(fofb) + abs(fofa)); // Linear Interpolation

        fofc = f(c);

        if(fofc < 0) {
            a = c;
        }
        else if(fofc == 0) {
            a = c;
            break;
        }
        else {
            b = c;
        }

    }

    printf("The root is %g, to %d decimal places, after %d iterations.\n", (double)a, dp, i);
}

Answer 1

函数 abs()（来自 <stdlib.h>）具有签名 int abs(int); — 您正在从计算中获取整数。

您应该使用 <math.h> 中的 long double fabsl(long double);。

您还应该使用 powl() 而不是 pow()（long double 与 double），并且 roundl() 而不是 roundf()（ long double 对比 float).

换句话说，请确保您使用的是正确的类型。

当您解决了类型问题后，您仍然会遇到收敛问题。如果您制作了 MCVE (Minimal, Complete, Verifiable Example)，那将会有所帮助。但是，这是我可以从您的问题中推断出的 MCVE：

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

static inline long double f(long double x)
{
    return(powl(x, 3) + 2 * x - 8);
}

int main(void)
{
    long double a = 1.0L;
    long double b = 2.0L;
    int dp = 6;

    while (roundl(a * powl(10, dp)) / powl(10, dp) != roundl(b * powl(10, dp)) / powl(10, dp))
    {
        long double fofa = f(a);
        long double fofb = f(b);
        long double c = (b * fabsl(fofa) + a * fabsl(fofb)) / (fabsl(fofb) + fabsl(fofa));
        long double fofc = f(c);
        printf("a = %+.10Lf, f(a) = %+.10Lf\n", a, fofa);
        printf("b = %+.10Lf, f(b) = %+.10Lf\n", b, fofb);
        printf("c = %+.10Lf, f(c) = %+.10Lf\n", c, fofc);
        putchar('\n');

        if (fofc < 0.0L)
        {
            a = c;
        }
        else if (fofc == 0.0L)
        {
            a = c;
            break;
        }
        else
        {
            b = c;
        }
    }
    printf("Result: a = %Lg\n", a);
    return 0;
}

我从中得到的输出是：

a = +1.0000000000, f(a) = -5.0000000000
b = +2.0000000000, f(b) = +4.0000000000
c = +1.5555555556, f(c) = -1.1248285322

a = +1.5555555556, f(a) = -1.1248285322
b = +2.0000000000, f(b) = +4.0000000000
c = +1.6531049251, f(c) = -0.1762579238

a = +1.6531049251, f(a) = -0.1762579238
b = +2.0000000000, f(b) = +4.0000000000
c = +1.6677455452, f(c) = -0.0258828049

a = +1.6677455452, f(a) = -0.0258828049
b = +2.0000000000, f(b) = +4.0000000000
c = +1.6698816424, f(c) = -0.0037639074

a = +1.6698816424, f(a) = -0.0037639074
b = +2.0000000000, f(b) = +4.0000000000
c = +1.6701919841, f(c) = -0.0005465735

a = +1.6701919841, f(a) = -0.0005465735
b = +2.0000000000, f(b) = +4.0000000000
c = +1.6702370440, f(c) = -0.0000793539

a = +1.6702370440, f(a) = -0.0000793539
b = +2.0000000000, f(b) = +4.0000000000
c = +1.6702435859, f(c) = -0.0000115206

a = +1.6702435859, f(a) = -0.0000115206
b = +2.0000000000, f(b) = +4.0000000000
c = +1.6702445357, f(c) = -0.0000016726

a = +1.6702445357, f(a) = -0.0000016726
b = +2.0000000000, f(b) = +4.0000000000
c = +1.6702446735, f(c) = -0.0000002428

a = +1.6702446735, f(a) = -0.0000002428
b = +2.0000000000, f(b) = +4.0000000000
c = +1.6702446936, f(c) = -0.0000000353

a = +1.6702446936, f(a) = -0.0000000353
b = +2.0000000000, f(b) = +4.0000000000
c = +1.6702446965, f(c) = -0.0000000051

a = +1.6702446965, f(a) = -0.0000000051
b = +2.0000000000, f(b) = +4.0000000000
c = +1.6702446969, f(c) = -0.0000000007

a = +1.6702446969, f(a) = -0.0000000007
b = +2.0000000000, f(b) = +4.0000000000
c = +1.6702446970, f(c) = -0.0000000001

a = +1.6702446970, f(a) = -0.0000000001
b = +2.0000000000, f(b) = +4.0000000000
c = +1.6702446970, f(c) = -0.0000000000

a = +1.6702446970, f(a) = -0.0000000000
b = +2.0000000000, f(b) = +4.0000000000
c = +1.6702446970, f(c) = -0.0000000000

a = +1.6702446970, f(a) = -0.0000000000
b = +2.0000000000, f(b) = +4.0000000000
c = +1.6702446970, f(c) = -0.0000000000

a = +1.6702446970, f(a) = -0.0000000000
b = +2.0000000000, f(b) = +4.0000000000
c = +1.6702446970, f(c) = -0.0000000000

无限循环的原因很清楚； a 和 b 之间的差异不是一个小分数。您需要检查您的收敛标准。它可能应该将 fofc 与 0.0 进行比较，以达到给定的小数位数 - 或者类似的东西。

Answer 2

您正在实施的方法称为 regula falsi 或 false position method。

只要保持相反的符号条件，实际上不需要使用绝对值。

有一个众所周知的 stalling 普通 vanilla regula falsi 问题，因为此时函数在剩余区间上是凸的，一个端点将不再移动向根。可以通过简单的修改来避免这种情况，例如插入对分步骤。 Illinois modification 实施起来更简单但更难理解。有关详细信息，请参阅维基百科关于正则假的文章。

或者这个问答：Regula-Falsi Algorithm?

改编自 link 中的答案：

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

long double f(long double x) {
   return powl(x, 3) + 2 * x - 8;
}

int main(void) {
    const int dp = 18;
    long double eps=0.5*powl(10,-dp);
    int i=0;
    long double a=1, fofa = f(a);
    long double b=2, fofb = f(b);

    printf("\na=%.21Lf b=%.21Lf  fofa=%.21Lf  fofb=%.21Lf\n------\n",a,b, fofa,fofb);

    if(signbit(fofb)==signbit(fofa)) {
        printf("Warning, initial values do not have opposite sign!\n");
    }
    do {
        long double c=(a*fofb-b*fofa)/(fofb-fofa), fofc = f(c);

        if( signbit(fofc)!=signbit(fofa) ) {
            b=a; fofb=fofa;
            a=c; fofa=fofc;
        } else {
            a=c; fofa=fofc;
            fofb *= 0.5;
        }
        i++;  
        printf("\na=c=%.21Lf b=%.21Lf  fofa=fofc=%.21Lf  fofb=%.21Lf",c,b, fofc,fofb);

    } while(fabsl(b-a)>eps);

    printf("\ngoal reached after %d iterations\n",i);

    return 0;
}

结果

a=1.000000000000000000000 b=2.000000000000000000000  fofa=-5.000000000000000000000  fofb=4.000000000000000000000
------

a=c=1.555555555555555555507 b=2.000000000000000000000  fofa=fofc=-1.124828532235939643549  fofb=2.000000000000000000000
a=c=1.715539947322212467064 b=1.555555555555555555507  fofa=fofc=0.480046589479470829469  fofb=-1.124828532235939643549
a=c=1.667685780603345490963 b=1.715539947322212467064  fofa=fofc=-0.026500999000164700194  fofb=0.480046589479470829469
a=c=1.670189362207942139265 b=1.715539947322212467064  fofa=fofc=-0.000573759143326624515  fofb=0.240023294739735414734
a=c=1.670297511133477909220 b=1.670189362207942139265  fofa=fofc=0.000547652143260468627  fofb=-0.000573759143326624515
a=c=1.670244695550498326532 b=1.670297511133477909220  fofa=fofc=-0.000000014643676336194  fofb=0.000547652143260468627
a=c=1.670244696962696986627 b=1.670297511133477909220  fofa=fofc=-0.000000000000373724489  fofb=0.000273826071630234313
a=c=1.670244696962769068724 b=1.670244696962696986627  fofa=fofc=0.000000000000373706274  fofb=-0.000000000000373724489
a=c=1.670244696962733028543 b=1.670244696962769068724  fofa=fofc=-0.000000000000000000434  fofb=0.000000000000373706274
a=c=1.670244696962733028651 b=1.670244696962733028543  fofa=fofc=0.000000000000000000867  fofb=-0.000000000000000000434
goal reached after 10 iterations

为什么 C 为这个算法赋予 Python 不同的值？

Why does C give different values to Python for this algorithm?

c

math

linear-interpolation