c++ - 使用 Runge-Kutta(二阶)积分波动方程

Question

我尝试求解简单的数值方程 - 没有源的线性波动方程:u_tt = v² u_xx

其中 v - 波速。

我使用初始条件:

u(x, 0) = sin(x)
u_x(x, 0) = -v * sin(x)

对应于沿 x 轴以速度 v 传播的初始波 sin(x)(sin(x-vt) 是该波动方程的解)。

因此，问题是在 2-3 次迭代期间，数值解与解析答案 (sin(x-vt)) 一致，但在第 3 次迭代后出现问题，解增长非常快(取值约为 128)。

我使用 Runge-Kutta 二阶排序法(欧拉法)和二阶导数的四阶近似。

这是一段 C++ 代码:

//2nd derivative function

float z2der(int n, int i, int j, float* x, float** zx)
{
    float h = x[2] - x[1];
    float zr, zl, zm, zrr, zll;

    zm = -30.0 * zx[i][j];
    if (i < n - 1)
        zr = 16.0 * zx[i + 1][j];
    else
        zr = 16.0;
    if (i > 0)
        zl = 16.0 * zx[i - 1][j];
    else
        zl = 16.0;
    if (i < n - 2)
        zrr = -zx[i + 2][j];
    else
        zrr = 1.0; //CHECK
    if (i > 1)
        zll = -zx[i - 2][j];

    if (i > 1 && i < n - 2)
        return (zll + zrr + zl + zr + zm) / (12.0 * h * h);
    if (i == 1)
        return (15 / 4 * zx[i][j] - 77 / 6 * zx[i + 1][j] - 107 / 6 * zx[i + 2][j] - 13 * zx[i + 3][j] - 61 / 12 * zx[i + 4][j]) / (h*h);
    if (i == n - 2)
        return (15 / 4 * zx[i][j] - 77 / 6 * zx[i - 1][j] - 107 / 6 * zx[i - 2][j] - 13 * zx[i - 3][j] - 61 / 12 * zx[i - 4][j]) / (h*h);
    if (i == 0 || i == n - 1)
        return 0;
}


//PARAMETERS

float tmin = 0.0;
float tmax = 12.0;
int nt = 37000;
float dt = (tmax - tmin) / float(nt);
float *time = new float[nt];
for (int j = 0; j < nt; j++)
    time[j] = tmin + dt * float(j);

float xmin = -30.0;
float xmax = 30.0;
int nx = 901;
float *x = new float[nx];
for (int i = 0; i < nx; i++)
    x[i] = xmin + (xmax - xmin) * float(i) / float(nx);

float v = -0.1; // velocity

float** wave = new float*[nx];
for (int i = 0; i < nx; i++)
    wave[i] = new float[nt];

float** waved = new float*[nx];
for (int i = 0; i < nx; i++)
    waved[i] = new float[nt];

for (int i = 0; i < nx; i++)
{
    xx = x[i];
    wave[i][0] = sin(xx);
    waved[i][0] = -v * v * sin(xx);
}

float* wavedd = new float[nx];
float* waveddot = new float[nx];

for (int i = 0; i < nx; i++) //j=0;
{
    wavedd[i] = z2der(nx, i, 0, x, wave);
    waveddot[i] = v * v * wavedd[i];
}

for (int j = 0; j < nt - 1; j++)
{
    for (int i = 0; i < nx; i++)
    {
        waved[i][j + 1] = waved[i][j] + dt * waveddot[i];
        wave[i][j + 1] = wave[i][j] + dt * waved[i][j];

        wavedd[i] = z2der(nx, i, j + 1, x, wave);
        waveddot[i] = v * v * wavedd[i];
    }
}

所以，我不知道为什么会出错，是什么原因。

Answer 1

你所有的 RK 画面值都会有问题，因为你正在进行整数除法

return (15 / 4 * ...

由于运算符的优先级，这将评估为

return (3 * ...

你需要完全在 double 中工作

return (15.0 / 4.0 * ...