数值解不会像它应该的那样发散。为什么？

Question

我正在尝试使用欧拉方法来逼近一个微分方程：

u'=3*(u-t)
u(0)=1/3

问题应该在具有浮点精度的大量步骤中发散。这是由于初始数据的舍入误差。 对我的一些朋友来说，这段代码有分歧，但对我来说却没有。

编译器会不会提高精度？

评论更新：

@KillzoneKid no g++ -Wall -pedantic main.cpp 不打印任何内容

@some-programmer-dude 实际输出是精确解，而由于浮动错误，输出应该发散。类似 7189 而不是 10+ 1/3（精确解）

@pac0 编译后的输出不起作用，因为其他人都在使用非常旧的 linux 版本（内核 2.6），但我会尝试设置编译器选项。

趋势是使用 mac 的人（我的朋友）和使用现代 linux 的人（我）会遇到这个问题，而使用 windows 或非常旧的 linux 的人不会。

@1201ProgramAlarm 我在 XPS 15 上使用 ubuntu 17.10，并使用 CLion 作为编辑器。 为简单起见，我使用的是操作系统中捆绑的 g++。

#include <iostream>
#include <math.h>
#include <fstream>

using namespace std;
typedef float Real;
Real f(Real t,Real u);
const Real pi=4.0*atan(1.0);

int main() {

    Real u,t,T,tau;


    // print to file
    char n_file[21]={0};
    std::cout << "file name: " << std::endl;
    std::cin >> n_file ;

    ofstream prt(n_file);
    prt.precision(15);
    //

    t=0.0;//start time
    T=10.0;//final time
    u=1.0/3.0;// u(0)
    unsigned long N;
    for(int k=1;k<=20;k++){
        N=(unsigned long)pow(10,k);
        tau=(T-t)/Real(N);
        for(unsigned long n=1;n<=N;n++){
            u+=tau*f(t,u);
            t+=tau;
        }
        prt << "With " << N << " steps, at time " << tau << "result is " <<u<<endl;

        prt << endl << endl << endl;

        u = 1.0/3.0;
        t = 0.0;

    }

//

    return 0;
}


//
Real f(Real t, Real u)
{
    return 3*(u-t);
}

Answer 1

u=1.0/3.0 + 1e-8; 的计算以双精度完成，然后在分配给 u 时四舍五入到最接近的浮点值。使用 Windows 和 Visual Studio，在从 double 到 float 的转换过程中，+1e-8 被舍入，但 +1e-7 大到足以影响 float 结果：

1.0/3.0 + 1e-8 == 1.0/3.0 == 32 bit hex integer 0x3eaaaaab
                          ~= 0.33333334

1.0/3.0 + 1e-7 ==            32 bit hex integer 0x3eaaaaae
                          ~= 0.33333343

---

环境之间的差异可能是浮点控制字中的舍入设置问题，也可能是硬件实现问题。

---

我将代码更改为使用 2 的幂（在本例中为 8 的幂）的步数，并将最大步数限制为对应于浮点数尾数部分的有效位数.这消除了分歧，部分原因是 u 的初始值略大于 1/3，而 t 和 tau 是精确的（因为步数是 2 的幂），部分原因是乘以 tau减少误差大于重复加法增加误差。将 u 初始化为 1/3 - 1e-7，总和在 64 步处发散，变为 -5770，但 8 步时为 10.30，>= 512 步时为 10.333333。

#include <iomanip>
#include <iostream>
#include <math.h>
#include <fstream>

using namespace std;
typedef float Real;
Real f(Real t,Real u);

int main() {
    Real u,t,T,tau;

    // print to file
    char n_file[21]={0};
    std::cout << "file name: " << std::endl;
    std::cin >> n_file ;

    ofstream prt(n_file);
    prt.precision(15);

    T=10.0;    //final time
    unsigned long N = 1;
    for(int k=1;k<=7;k++){
        N *= 8;                 // # steps is power of 2
        u = (Real)(1.0/3.0);
        t = 0.0;
        tau=T/Real(N);
        for(unsigned long n=1;n<=N;n++){
            u+=tau*f(t,u);
            t+=tau;
        }
        prt << "With " << setw(8) << N << " steps result is " << u <<endl;
    }

    return 0;
}

Real f(Real t, Real u)
{
    return 3*(u-t);
}