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

Question

我正在尝试使用欧拉方法来近似微分方程:

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

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

编译器是否有可能提高精度？

评论更新:

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

@some-programmer-dude 实际输出是精确的解决方案，而由于 float 错误，输出应该有所不同。像 7189 而不是 10+ 1/3(精确解)

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

趋势是使用 mac 的人(我的 friend )和使用现代 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; 的计算以 double 完成，然后在分配给 u 时四舍五入为最接近的浮点值。使用 Windows 和 Visual Studio，+1e-8 在从 double 到 float 的转换过程中四舍五入，但 +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 的幂)的步数，并将最大步数限制为与 float 的尾数部分中的有效位数相对应.这消除了分歧，部分原因是 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);
}