Regula-Falsi算法是什么？

Question

Regula-Falsi算法是什么？

4

我正在尝试实现Regula-Falsi算法来解决方程2(x^3)-x-2的问题，但问题是变量c的值保持不变，即使我的代码应该改变它。

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


float fonc(float x)
{
    int result;
    result=2*(pow(x,3))-x-2;
    return result;
}

int main(void)
{
    float eps=pow(10,-4);
    int i=0;
    float a,b,c;
    a=1;
    b=2;
    do
    {
        c=((a*fonc(b))-(b*fonc(a)))/((fonc(b)-fonc(a)));
        if(fonc(c)*fonc(a)<0)
        {
            b=c;
        }
        else
        {
            a=c;
        }
        i++;    
        printf("\n%f",c);
    }
    while(fabs(b-c)>eps);

    printf("le nombre d'itération %d",i);
}

- user3396807

建议使用 double 而不是 float。 - pmg

即使使用双精度，我们也会到达某个无法正确获取精度的地方。 - user3396807

我的意思是一般而言，不是特指这个程序。如果你需要一个浮点数，你的第一反应应该是使用 double。 - pmg

不要使用 pow(10,-4)，请改用 1e-4 或 1e-4f。 - phuclv

3个回答

3

即使所有数据类型都正确，按照所呈现的算法可能出现哪些问题?

与二分法不同，纯正则法不会强制使区间长度趋近于零。如果使用单调且凸函数，则迭代将在仅改变区间一侧的情况下停止，并呈几何收敛。任何充分平滑的函数最终都会在剩余的括号内具有这些属性。

为了正确捕捉这种行为，在计算后应立即将中点与区间端点,进行比较。作为奖励分数，检查处的值是否足够小，如果为真，无论离区间的端点有多远，也应该中断迭代。

存在许多简单的技巧来强制使区间长度为零。复杂的技巧导致Brent方法。其中一个简单的技巧是伊利诺伊变体。在这些变体中，中点被认为是一个凸和。

c = |f(b)|/(|f(a)|+|f(b)|) * a + |f(a)|/(|f(a)|+|f(b)|) * b

由于 f(a) 和 f(b) 的符号相反，因此这等效于原始公式。如果侧边 b 不变，则在减小函数值 f(b) 的情况下增加其在这个凸和中的重要性，即将其乘以附加权重因子。这将使中点 c 向 b 移动，在很少的步骤中找到一个将替换 b 的中点。
以下是实现 regula falsi（或者说是 false position 方法）的 Illinois 变体。该算法在 6 次迭代中找到了函数值为 2.2e-6、包含区间长度为 6e-7 的解。
#include<math.h> #include<stdio.h> float fonc(float x) { return (2*x*x-1)*x-2; } int main(void) { float eps=1e-6; int i=0; float a=1, fa = fonc(a); float b=2, fb = fonc(b); printf("\na=%10.7f b=%10.7f fa=%10.7f fb=%10.7f\n------\n",a,b, fa,fb); if(signbit(fb)==signbit(fa)) { printf("Attention, les valeurs initiales de 'fonc' n'ont pas de signe opposeés!\n"); } do { float c=(a*fb-b*fa)/(fb-fa), fc = fonc(c); if( signbit(fc)!=signbit(fa) ) { b=a; fb=fa; } else { fb *= 0.5; } a=c; fa=fc; i++; printf("\na=c=%10.7f b=%10.7f fa=fc=%10.7f fb=%10.7f",c,b, fc,fb); if(fabs(fc)<eps) break; } while(fabs(b-a)>eps); printf("\nle nombre d'itération %d\n",i); return 0; }

输出结果为
a= 1.0000000 b= 2.0000000 fa=-1.0000000 fb=12.0000000 ------ a=c= 1.0769231 b= 2.0000000 fa=fc=-0.5789710 fb= 6.0000000 a=c= 1.1581569 b= 2.0000000 fa=fc=-0.0512219 fb= 3.0000000 a=c= 1.1722891 b= 1.1581569 fa=fc= 0.0497752 fb=-0.0512219 a=c= 1.1653242 b= 1.1722891 fa=fc=-0.0003491 fb= 0.0497752 a=c= 1.1653727 b= 1.1722891 fa=fc=-0.0000022 fb= 0.0248876 a=c= 1.1653733 b= 1.1653727 fa=fc= 0.0000020 fb=-0.0000022 le nombre d'itération 6

- Lutz Lehmann

2

一个问题是 result 是 int 类型。你几乎肯定希望它是 float 或者 double，否则 fonc() 的结果会被截断为整数。

- NPE

尝试更改，问题仍然存在。 - user3396807

@AbdelhadiKhiati：很奇怪，这应该可以解决它。你保存了文件吗？重新编译/链接了吗？运行了正确的二进制文件吗？ - NPE

@AbdelhadiKhiati 这真的解决了所有问题吗？在我的端上，仅有这个更改仍然会导致程序进入无限循环。当然，这个更改是绝对必要的，但代码似乎仍需要更多关注。 - Utkan Gezer

@ThoAppelsin 是的，这对我做了一些改变，循环变成了无限循环，有时则是有限循环。我真的不知道它在发生什么，但无论我如何改变(eps)的精度，结果都是稳定的。:/ - user3396807

@AbdelhadiKhiati 参考我刚才给出的答案。另外，为了美观起见，您可能需要将 printf 中的 "\n%f" 更改为 "%f\n"。 - Utkan Gezer

如果您写出迭代的更多值（例如，打印出所有a，fonc(a)，b，fonc(b)），您将会对实现有更深入的了解。为了保持输出简洁，在循环中插入if(i>40) break;和if(fabs(fonc(c))<eps) break;，作为第一条或最后一条语句。 - Lutz Lehmann

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- Utkan Gezer · Accepted Answer

好的，为了纠正算法以实现其本应实现的功能：

除了将result的类型从int改为float外，您还需要更改循环条件。目前是这样的：

    while ( fabs( b - c ) > eps );

这意味着循环将继续进行，直到值b和c之间的距离低于0.0001，但在此条件下，当前代码至少在我的设备上会一直运行。

无论如何，我们并不想减少b和c之间的差异。我们真正想要的是使fonc(c) = 2*c*c*c - c - 2小于我们的eps。最终，我们希望它尽可能接近零，以便c大致成为该函数的一个根。因此简单地说：

    while ( fabs( fonc( c ) ) > eps );

这就是我们应该的情况。这样，连同 int --> float 的更改，它不会陷入无限循环，在14次迭代中完成工作。