Численное решение уравнения синус-Гордона

[lv00]

Здравствуйте.

Я пытаюсь решить численно одномерное уравнение синус-Гордона

$$\phi_{tt} - \phi_{xx} + \sin \phi = 0 \ \, $$

а именно, воспроизвести эволюцию движущегося кинка

$$\phi_k (x, t) = 4 \ arctan \ e^{(x-vt)/ \sqrt{1-v^2}} $$

Для решения использую сеточную аппроксимацию $\phi_{xx}$ и явный метод Рунге-Кутты 4-го порядка точности для получения временной эволюции.
Но проблема состоит в том, что кинк распадается.

На этом рисунке показан кинк в начальный момент времени $t=0$.

На этом показано значение функции $\phi(x)$ при $t=20$.


Привожу код:

Используется синтаксис C++

//Аппроксимация второй производной//

double z2der(int n, int i, int j, double* x, double** zx) {
        double h = x[2] - x[1];
        double zr, zl, zm, zrr, zll;
    
        if (i > 1 && i<n - 2) {
                zr = 16.0*zx[i + 1][j];
                zl = 16.0*zx[i - 1][j];
                zrr = -zx[i + 2][j];
                zll = -zx[i - 2][j];
        zm = -30.0*zx[i][j];
                return (zll + zrr + zl + zr + zm) / (12.0 * h * h);
    }
    if (i==0 || i==n-1)
        return 0;
    if (i==1 || i==n-2)
        return 0;
}

double f_rhs(int i, int j, double **arr) {
    return sin(arr[i][j]); // Sine-Gordon
}

int main() {
//Параметры решетки//

        double xmin = -60.0;
        double xmax = 60.0;
        int nx = 1000;
        double *x = new double[nx];
        for (int i = 0; i<nx; i++)
                x[i] = xmin + (xmax - xmin)*double(i) / double(nx);
    
        double tmin = 0.0;
        double tmax = 100.0;

int nt = 10000;
        double dt = (tmax - tmin) / double(nt);
        double *time = new double[nt];
        for (int j = 0; j<nt; j++)
                time[j] = tmin + dt*double(j);

double v = 0.30; // скорость кинка
double a = 10; //начальный сдвиг кинка

//4th order RK method
    
        double **k1 = new double*[nx]; //for phi
        for (int i = 0; i < nx; i++)
                k1[i] = new double[nt];
        double **k2 = new double*[nx];
        for (int i = 0; i < nx; i++)
                k2[i] = new double[nt];
        double **k3 = new double*[nx];
        for (int i = 0; i < nx; i++)
                k3[i] = new double[nt];
        double **k4 = new double*[nx];
        for (int i = 0; i < nx; i++)
                k4[i] = new double[nt];
    
    
        double **q1 = new double*[nx]; //for phi_t
        for (int i = 0; i < nx; i++)
                q1[i] = new double[nt];
        double **q2 = new double*[nx];
        for (int i = 0; i < nx; i++)
                q2[i] = new double[nt];
        double **q3 = new double*[nx];
        for (int i = 0; i < nx; i++)
                q3[i] = new double[nt];
        double **q4 = new double*[nx];
        for (int i = 0; i < nx; i++)
                q4[i] = new double[nt];
    
        double **r2 = new double*[nx];
        for (int i = 0; i < nx; i++)
                r2[i] = new double[nt];
        double **r3 = new double*[nx];
        for (int i = 0; i < nx; i++)
                r3[i] = new double[nt];
        double **r4 = new double*[nx];
        for (int i = 0; i < nx; i++)
                r4[i] = new double[nt];
    
        double** phi = new double*[nx];
        for (int i = 0; i < nx; i++)
                phi[i] = new double[nt];
    
        double** phi_t = new double*[nx];
        for (int i = 0; i < nx; i++)
                phi_t[i] = new double[nt];

double xv1; //вспомогательная координата (см. далее)

//Начальные условия//

for (int i = 0; i < nx; i++) {
                xx = x[i];
                xx1 = x[i] + a;
        xv1 = xx1/sqrt(1.0-v*v);

 // --- SINE-GORDON --- //
        phi[i][0] = 4*atan(exp(xv1)); //сам кинк
        phi_t[i][0] = (-v/sqrt(1-v*v))*4.0*exp(xv1)/(1+pow(exp(xv1), 2.0)); //его производная по времени
}


//Решение//
for (int j = 0; j < nt; j++) {
                tt = j * dt;

for (int i = 0; i < nx; i++) {
                        k1[i][j] = phi_t[i][j];
                        q1[i][j] = (z2der(nx, i, j, x, phi) + f_rhs(i, j, phi));
                }
                for (int i = 0; i < nx; i++) {
                        r2[i][j] = phi[i][j] + dt * k1[i][j] / 2.0;
                        k2[i][j] = phi_t[i][j] + dt * q1[i][j] / 2.0;
                }
                for (int i = 0; i < nx; i++) {
                        q2[i][j] = (z2der(nx, i, j, x, r2) + f_rhs(i, j, r2));
                        r3[i][j] = phi[i][j] + dt * k2[i][j] / 2.0;
                        k3[i][j] = phi_t[i][j] + dt * q2[i][j] / 2.0;
                }
                for (int i = 0; i < nx; i++) {
                        q3[i][j] = (z2der(nx, i, j, x, r3) + f_rhs(i, j, r3));
                        r4[i][j] = phi[i][j] + dt * k3[i][j];
                        k4[i][j] = phi_t[i][j] + dt * q3[i][j];
                }
                for (int i = 0; i < nx; i++)
                        q4[i][j] = (z2der(nx, i, j, x, r4) + f_rhs(i, j, r4));
        for (int i = 1; i < nx - 1; i++) {
                        phi[i][j + 1] = phi[i][j] + dt / 6.0 * (k1[i][j] + 2.0 * k2[i][j] + 2.0 * k3[i][j] + k4[i][j]);
                        phi_t[i][j + 1] = phi_t[i][j] + dt / 6.0 * (q1[i][j] + 2.0 * q2[i][j] + 2.0 * q3[i][j] + q4[i][j]);
                }
        }
 }
 

[TelmanStud]

Вот в этой небольшой брошюрке
http://pauli.uni-muenster.de/tp/fileadmin/lehre/NumMethoden/WS1011/script1011.pdf
использована самая простая явная схема, и вроде кинк "нормально" бежит

[lv00]

TelmanStud в сообщении #1060770 писал(а):

Вот в этой небольшой брошюрке
http://pauli.uni-muenster.de/tp/fileadmin/lehre/NumMethoden/WS1011/script1011.pdf
использована самая простая явная схема, и вроде кинк "нормально" бежит

Это замечательно, только там метод основан на решении систем линейных уравнений. Насколько я знаю по опыту, обычно с такими вещами проблем не возникает.
В моем же случае, по всей видимости, неправильно работает сам метод Рунге-Кутты, но в чем дело - я не могу понять.

[lv00]

UPD: заметил, что не скопировал фиксированные граничные условия для кинка: они в разделе SOLUTION задаются вручную следующим образом:

Используется синтаксис C++

 
                phi[0][j] = 0.0; 
                phi[nx - 1][j] = 6.28319;
 

[TelmanStud]

lv00 в сообщении #1060858 писал(а):

Это замечательно, только там метод основан на решении систем линейных уравнений.

Там же вроде явная схема, какие системы уравнений?

[lv00]

TelmanStud в сообщении #1061984 писал(а):

lv00 в сообщении #1060858 писал(а):

Это замечательно, только там метод основан на решении систем линейных уравнений.

Там же вроде явная схема, какие системы уравнений?

Я имел в виду, что там используются не интеграторы типа Рунге-Кутты, а схемы конечных разностей (могут быть явными и неявными) - в общем случае нужно решать систему линейных уравнений.

[Pphantom]

Так вообще-то явные схемы второго и более высоких порядков монотонными не бывают. Не в этом ли дело?