2014 dxdy logo

Научный форум dxdy

Математика, Физика, Computer Science, Machine Learning, LaTeX, Механика и Техника, Химия,
Биология и Медицина, Экономика и Финансовая Математика, Гуманитарные науки




Начать новую тему Ответить на тему
 
 Численное решение уравнения синус-Гордона
Сообщение09.10.2015, 13:04 


18/05/14
71
Здравствуйте.

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

$$\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]);
                }
        }
 }
 

 Профиль  
                  
 
 Re: Численное решение уравнения синус-Гордона
Сообщение09.10.2015, 14:05 
Аватара пользователя


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

 Профиль  
                  
 
 Re: Численное решение уравнения синус-Гордона
Сообщение09.10.2015, 18:24 


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

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

 Профиль  
                  
 
 Re: Численное решение уравнения синус-Гордона
Сообщение12.10.2015, 14:49 


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

Используется синтаксис C++
 
                phi[0][j] = 0.0;
                phi[nx - 1][j] = 6.28319;
 

 Профиль  
                  
 
 Re: Численное решение уравнения синус-Гордона
Сообщение13.10.2015, 13:38 
Аватара пользователя


05/04/13
580
lv00 в сообщении #1060858 писал(а):
Это замечательно, только там метод основан на решении систем линейных уравнений.

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

 Профиль  
                  
 
 Re: Численное решение уравнения синус-Гордона
Сообщение13.10.2015, 16:19 


18/05/14
71
TelmanStud в сообщении #1061984 писал(а):
lv00 в сообщении #1060858 писал(а):
Это замечательно, только там метод основан на решении систем линейных уравнений.

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


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

 Профиль  
                  
 
 Re: Численное решение уравнения синус-Гордона
Сообщение13.10.2015, 17:59 
Заслуженный участник


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

 Профиль  
                  
Показать сообщения за:  Поле сортировки  
Начать новую тему Ответить на тему  [ Сообщений: 7 ] 

Модераторы: Karan, Toucan, PAV, maxal, Супермодераторы



Кто сейчас на конференции

Сейчас этот форум просматривают: нет зарегистрированных пользователей


Вы не можете начинать темы
Вы не можете отвечать на сообщения
Вы не можете редактировать свои сообщения
Вы не можете удалять свои сообщения
Вы не можете добавлять вложения

Найти:
Powered by phpBB © 2000, 2002, 2005, 2007 phpBB Group