Интегрирование волнового уравнения с Рунге-Кутта (2-й порядок)

Я пытаюсь решить численно простое уравнение — линейное волновое уравнение без источников: u_тт = v² U_хх

где v — скорость волны.

Я использую начальные условия:

и (х, 0) = грех (х)
U_Икс(х, 0) = -v * грех (х)

которые соответствуют начальной волне sin (x), распространяющейся со скоростью v вдоль оси x (sin (x-vt) является решением этого волнового уравнения).

Таким образом, проблема состоит в том, что в течение 2-3 итераций численное решение совпадает с аналитическим ответом (sin (x-vt)), но после 3-й итерации что-то идет не так, и решение очень быстро растет (принимает значения около 128).

Я использую метод Рунге-Кутты 2-го рудера (метод Эйлера) с приближением 4-го порядка 2-й производной.

Вот код на 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];
}
}

Итак, я не знаю, почему это идет не так и какова причина.