数値的に単純な方程式を解こうとしています-ソースのない線形波動方程式: u tt = v 2 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 の値を取ります)。
ルンゲ・クッタの二次導関数法（オイラー法）を用いて二次導関数を四次近似します。
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];
}
}
だから、なぜそれがうまくいかないのか、その理由は何なのかわかりません。