Regula-Falsi 算法？

Question

我正在尝试实现 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);
}

Answer 1

好的，只是为了更正算法以实际执行它应该执行的操作：

除了将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 次迭代中完成工作。

Answer 2

添加了前言：即使所有数据类型都合适，所呈现的算法也会出现什么问题？

与二分法相比，纯regula falsi 不会强制区间长度为零。如果使用单调和凸函数，迭代将停滞在仅改变间隔的一侧并具有几何收敛性。

为了正确捕捉这种行为，中点 c 应在计算后立即与两个区间末端进行比较。作为奖励点，检查 c 处的值是否足够小，如果为真，也中断迭代，不管它离区间末端有多远。

有许多简单的技巧可以强制间隔长度为零。复杂的技巧导致了布伦特的方法。一个简单的技巧是伊利诺伊变体。在这些变体中，中点被视为凸和

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

因为a f(a) 和f(b) 的符号相反，所以这等价于原来的公式。如果边 b 没有变化，则通过减小函数值 f(b)，即乘以额外的权重因子，增加其在该凸和中的重要性。

---

以下是 regula falsi（或错误位置方法）的伊利诺伊变体的实现。该算法在 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;
            a=c; fa=fc;
        }
        else
        {
            a=c; fa=fc;
            fb *= 0.5;
        }
        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

Answer 3

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