Решение смешанной задачи для уравнения синус-Гордона
09.01.2013, 15:40. Показов 2066. Ответов 3
Здравствуйте!
Требуется численно решить смешанную задачу для уравнения синус-Гордона, заданного на определённом отрезке и временном промежутке:
с начальными и краевыми условиями, благодаря которым задача имеет солитонное решение
Разностные схемы (явная и неявные) составлены и реализованы на C++:
явная
Кликните здесь для просмотра всего текста
| C++ | 1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
| #include <iostream>
#include <vector>
#define _USE_MATH_DEFINES
#include <cmath>
#include <fstream>
using namespace std;
#define C_L 10. // COORD_LENGTH - половина длины отрезка, на которой рассматривается уравнение ([-C_L, C_L])
#define C_S 0.1 // COORD_SPLIT - шаг по пространственной координате
#define T_L 20. // TIME_LENGTH - промежуток времени, на котором рассматривается уравнение ([0, T_L])
#define T_S 0.01 // TIME_SPLIT - шаг по временной координате
#define SPEED 1. // переменная уравнения - скорость волны
#define OMEGA 1. // переменная уравнения
// Вспомогательный макрос - запись вектора в файл
#define vout(u) for(int i = 0; i != (u).size(); i++)\
{\
input_file << (u)[i] << " ";\
}\
input_file << endl
int main()
{
vector<double> u_prev(2 * C_L / C_S), // следующая итерация
u_cur(2 * C_L / C_S), // текущая итерация
u_next(2 * C_L / C_S); // предыдущая итерация
ofstream input_file("D:\\simple_obvious.txt"); // файловый поток
// Установка начальных условий
for(int i = 0; i != u_prev.size(); i++)
{
u_prev[i] = 4 * atan(exp(-C_L + i * C_S));
u_cur[i] = u_prev[i] - T_S * 4 * exp(-C_L + i * C_S) / (1 + exp(2 * (-C_L + i * C_S)));
}
for(int j = 2; j <= (int)(T_L / T_S); j++)
{
// Установка краевых условий
u_next[0] = 4 * atan(exp(-C_L - SPEED * j * T_S));
u_next[u_next.size() - 1] = 4 * atan(exp(C_L - SPEED * j * T_S));
// Явная разностная схема
for(int i = 1; i != u_next.size() - 1; i++)
{
u_next[i] = 2 * u_cur[i] - u_prev[i] + (T_S * T_S) * (SPEED * SPEED) * (u_cur[i + 1] - 2 * u_cur[i] + u_cur[i - 1]) / (C_S * C_S) - (T_S * T_S) * (OMEGA * OMEGA) * sin(u_cur[i]);
}
// Запись в файл
vout(u_prev);
for(int i = 0; i != u_prev.size(); i++)
{
u_prev[i] = u_cur[i];
u_cur[i] = u_next[i];
u_next[i] = 0;
}
}
// Запись в файл
vout(u_cur);
return 0;
} |
|
и неявная
Кликните здесь для просмотра всего текста
| C++ | 1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
| #include <iostream>
#include <vector>
#define _USE_MATH_DEFINES
#include <cmath>
#include <fstream>
using namespace std;
#define C_L 10. // COORD_LENGTH - половина длины отрезка, на которой рассматривается уравнение ([-C_L, C_L])
#define C_S 0.1 // COORD_SPLIT - шаг по пространственной координате
#define T_L 20. // TIME_LENGTH - промежуток времени, на котором рассматривается уравнение ([0, T_L])
#define T_S 0.01 // TIME_SPLIT - шаг по временной координате
#define SPEED 1. // переменная уравнения - скорость волны
#define OMEGA 1. // переменная уравнения
// Вспомогательный макрос - запись вектора в файл
#define vout(u) for(int i = 0; i != (u).size(); i++)\
{\
input_file << (u)[i] << " ";\
}\
input_file << endl
// Процедура с методом прогонки (код взят из Википедии)
void procedure(vector<double> &, vector<double> &, vector<double> &);
int main()
{
vector<double> u_prev(2 * C_L / C_S), // следующая итерация
u_cur(2 * C_L / C_S), // текущая итерация
u_next(2 * C_L / C_S); // предыдущая итерация
ofstream input_file("D:\\simple_unobvious.txt"); // файловый поток
// Установка начальных условий
for(int i = 0; i != u_prev.size(); i++)
{
u_prev[i] = 4 * atan(exp(-C_L + i * C_S));
u_cur[i] = u_prev[i] - T_S * 4 * exp(-C_L + i * C_S) / (1 + exp(2 * (-C_L + i * C_S)));
}
for(int j = 2; j <= (int)(T_L / T_S); j++)
{
// Установка краевых условий
u_next[0] = 4 * atan(exp(-C_L - SPEED * j * T_S));
u_next[u_next.size() - 1] = 4 * atan(exp(C_L - SPEED * j * T_S));
// Применение метода прогонки
procedure(u_next, u_cur, u_prev);
// Запись в файл
vout(u_prev);
// Перенос данных для вычисления следующей итерации
for(int i = 0; i != u_prev.size(); i++)
{
u_prev[i] = u_cur[i];
u_cur[i] = u_next[i];
u_next[i] = 0;
}
}
// Запись в файл
vout(u_cur);
return 0;
}
void procedure(vector<double> &u_next, vector<double> &u_cur, vector<double> &u_prev)
{
int n = u_next.size() - 2; // число неизвестных
double m; // вспомогательная переменная
vector<double> a(n), // диагональ, лежащая под главной (нумеруется: [1;n-1])
b(n), // диагональ, лежащая над главной (нумеруется: [0;n-2])
c(n), // главная диагональ матрицы (нумеруется: [0;n-1])
f(n), // правая часть (столбец)
x(n); // решение, вектор x будет содержать ответ
// Установка необходимых значений
for(int i = 0; i < n; i++)
{
if(i > 0)
{
a[i] = 1.;
}
else
{
a[i] = 0.;
}
if(i < n - 1)
{
b[i] = 1.;
}
else
{
b[i] = 0.;
}
/* // Схема, полученная при sigma = 1/2
c[i] = -2 - 2 * (C_S * C_S) / (SPEED * SPEED * T_S * T_S);
f[i] = - 4. * (C_S * C_S) / (SPEED * SPEED * T_S * T_S) * u_cur[i + 1]
+ 2. * (C_S * C_S) / (SPEED * SPEED * T_S * T_S) * u_prev[i + 1]
- u_prev[i + 2] + 2 * u_prev[i + 1] - u_prev[i]
+ 2. * (OMEGA * OMEGA) * (C_S * C_S) / (SPEED * SPEED) * sin(u_cur[i + 1]);
*/
// Схема, полученная при sigma = 1/4
c[i] = -2. - 4. * (C_S * C_S) / (SPEED * SPEED * T_S * T_S);
f[i] = -2. * (u_cur[i + 2] - 2 * u_cur[i +1] + u_cur[i])
- u_prev[i + 2] + 2 * u_prev[i + 1] - u_prev[i]
+ 4. * (C_S * C_S) / (SPEED * SPEED * T_S * T_S) * (u_prev[i + 1] - 2 * u_cur[i + 1])
+ 4. * (C_S * C_S ) * (OMEGA * OMEGA) / (SPEED * SPEED) * sin(u_cur[i + 1]);
}
f[0] -= u_next[0];
f[n - 1] -= u_next[n + 1];
// Метод прогонки из Википедии
for (int i = 1; i < n; i++)
{
m = a[i] / c[i - 1];
c[i] = c[i] - m * b[i - 1];
f[i] = f[i] - m * f[i - 1];
}
x[n - 1] = f[n - 1] / c[n - 1];
for (int i = n - 2; i >= 0; i--)
{
x[i]=(f[i] - b[i] * x[i + 1]) / c[i];
}
// Запись данных в наш вектор
for(int i = 0; i < n; i++)
{
u_next[i + 1] = x[i];
}
return;
} |
|
.
Скрипты на python:
для явной схемы
Кликните здесь для просмотра всего текста
| Python | 1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
| #!/usr/bin/python
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation
with open("D:\\simple_obvious.txt") as f:
data = f.read()
data = data.split(' \n')
y = [row.split(' ') for row in data]
x = []
#z = []
#w = []
for i in range(len(y[0])):
x.append(i)
#for i in range(len(y) - 1):
# z.append(i)
# w.append(y[i][len(y[i]) - 1])
fig = plt.figure()
ax = fig.add_subplot(111)
line, = ax.plot(x, y[0])
ax.set_title("Plot title...")
ax.set_xlabel('your x label..')
ax.set_ylabel('your y label...')
# рисуем белым цветом лишь для того, чтобы анимация не выходила за пределы графика
for i in range(len(y) - 1):
ax.plot(x, y[i], c=str(1.))
# ax.plot(z, w, c = 'b')
def animate(i):
line.set_ydata(y[i]) # update the data
return line,
def init():
line.set_ydata(np.ma.array(y[0]))
return line,
ani = animation.FuncAnimation(fig, animate, np.arange(1, len(y) - 1), init_func=init,
interval=1, blit=True)
plt.show() |
|
и неявной
Кликните здесь для просмотра всего текста
| Python | 1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
| #!/usr/bin/python
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation
with open("D:\\simple_unobvious.txt") as f:
data = f.read()
data = data.split(' \n')
y = [row.split(' ') for row in data]
x = []
#z = []
#w = []
#v = []
for i in range(len(y[0])):
x.append(i)
#for i in range(len(y) - 1):
# z.append(i)
# w.append(y[i][len(y[0]) - 1])
# v.append(y[i][len(y[0]) - 2])
fig = plt.figure()
ax = fig.add_subplot(111)
line, = ax.plot(x, y[0])
ax.set_title("Plot title...")
ax.set_xlabel('your x label..')
ax.set_ylabel('your y label...')
#ax.plot(z, w, c='b')
#ax.plot(z, v, c='r')
for i in range(len(y) - 1):
ax.plot(x, y[i], c=str(1.))
def animate(i):
line.set_ydata(y[i]) # update the data
return line,
def init():
line.set_ydata(np.ma.array(y[0]))
return line,
ani = animation.FuncAnimation(fig, animate, np.arange(1, len(y) - 1), init_func=init,
interval=1, blit=True)
plt.show() |
|
рисуют итерации (анимация).
Получаем почти что идентичные изображения (отличия в итерациях ≈ 1e-4), но они не верны.
Не понимаю, почему итерации не сходятся к солитонному решению.
0
|
(
Сообщение от Genghis
