Решение краевой задачи

@dux99 · Регистрация: 03.05.2016

Студворк — интернет-сервис помощи студентам

Добрый вечер, уважаемые форумчане! Прошу помощи с решением проблемы: необходимо решить волновое уравнение с помощью явной разнотстной схемы типа "крестик":

$https://www.cyberforum.ru/cgi-bin/latex.cgi?\frac{\partial^2u(x,t)}{\partial t^{2}} = a^{2}\frac{\partial^2u(x,t)}{\partial x^{2}}$

при следующих начальных и граиничных условиях:

$https://www.cyberforum.ru/cgi-bin/latex.cgi?U(x,0)=0.2(1-x)sin(\pi x) \\ \frac{\partial U(x,0)}{\partial t} = 0 \\U(0,t) = 0 \\U(1,t) = 0$

Я попытался написать программу, которая позволяет изменять как граничные условия, делать анимированный и не анимированный график решения, сохранять значения расчетов в текстовый файл, но у меня возникает ошибка, до этого еще было куча ошибок, но как исправить и получить работающую программу я не знаю, я понимаю смысл ошибки, но не могу понять, почему такое решение здесь не применимо. Привожу пример кода на Python 3.9 и описание самой ошибки....

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
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
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
# Author: Alm99-collab
# Date: 25/10/2021 2:00
# Description: MSUT STANKIN Μyagkov Alexandr IDM - 21 - 03
 
 
"""
Программа для численного решения волнового уравнения в одномерном случае,
при разлинчых начальных и граничных условий, с помощью явной разностной
схемы типа "крест".
 
Общий вид волнового уравнения:
u_tt = a**2*u_xx + f(x,t) (0,L) где u=0 на диапазоне x=0,L, for t in (0,T].
Начальные условия в общем случае: u=I(x), u_t=V(x).
В случае неоднородного уравнения задается функция f(x,t).
 
Синтакис основной функции решателя:
u, x, t = sol(I, V, f, a, L, dt, C, T, user_func) где:
I = I(x) - функция.
V = V(x) - функция.
f = f(x,t) - функция.
U_0, U_L, - условия на границе.
C - число Куранта (=c*dt/dx), зависящее от шага dx. Является криетрием
стабильности численных расчето, если соблюдено условие (<=1)
dt - шаг по времени.
dx - шаг по координате.
T - конечное время симуляции волнового процесса.
user_func - функция (u, x, t, n) для вызова пользовательских сценариев,
таких как анимация или вывод графика, запси данных в текстовый файл,
расчет ошибки (в случае если известно точное решение) и.т.д.
"""
 
import numpy as np
import math
import matplotlib.pyplot as plt
import os
import time
import glob
 
 
def sol(I, V, f, a, L, C, T, U_0, U_L, dt, user_func = None):
    """
    Функция решатель волнового 1-D волнового уравнения
    u_tt = a**2*u_xx + f(x,t) (0,L) где u=0 на диапазоне
    x=0,L, for t in (0,T].
    :param I:
    :param V:
    :param f:
    :param a:
    :param L:
    :param C:
    :param T:
    :param U_0:
    :param U_L:
    :param dt:
    :param user_func:
    :return:
    """
 
    nt = int(round(T / dt))
    t = np.linspace(0, nt * dt, nt + 1)  # Массив точек для сетки времени
    dx = dt * a / float(C)
    nx = int(round(L / dx))
    x = np.linspace(0, L, nx + 1)  # Массив точек для сетки координат
    q = a ** 2
    C2 = (dt / dx) ** 2
    dt2 = dt * dt
 
    #  --- Проверка размерностей массивов сетки и шага между двумя точками сетки ---
    dx = x[1] - x[0]
    dt = t[1] - t[0]
 
    # --- Условие для проверки в качестве входного аргумента однородной функии ---
    if f is None or f == 0:
        f = lambda x, t: 0  # без использования векторизации используем лямбда- функцию
 
    # --- Условие для проверки в качестве входного аргумента начального условия dU(x,0)/dt ---
    if V is None or V == 0:
        V = lambda x: 0
 
    # Граничные условия для заданного волнового уравнения
    if U_0 is not None:
        if isinstance(U_0, (float, int)) and U_0 == 0:
            U_0 = lambda t: 0
    if U_L is not None:
        if isinstance(U_L, (float, int)) and U_L == 0:
            U_L = lambda t: 0
 
    # ---  Выделяем память под значения решений  ---
    u = np.zeros(nx + 1)  # Массив решений в узлах сетки на  временном шаге u(i,n+1)
    u_n = np.zeros(nx + 1)  # Массив решений в узлах сетки на  временном шаге u(i,n)
    u_nm = np.zeros(nx + 1)  # Массив решений в узлах сетки на  временном шаге u(i,n-1)
 
    # --- Проверка индексов для соблюдения размерностей массивов ---
    Ix = range(0, Nx + 1)
    It = range(0, Nt + 1)
 
    # --- Запись начальных условий ---
    for i in range(0, Nx + 1):
        u_n[i] = I(x[i])
 
    if user_func is not None:
        user_action(u_n, x, t, 0)
 
    # --- Разностная формулма явной схемы "типа крест" на первом шаге ---
    for i in Ix[1:-1]:
        u[i] = u_n[i] + dt * V(x[i]) + 0.5 * C2 * (0.5 * (q[i] + q[i + 1]) * (u_n[i + 1] - u_n[i]) -
                                                   0.5 * (q[i] + q[i - 1]) * (u_n[i] - u_n[i - 1])) + 0.5 * dt2 * f(
                x[i], t[0])
 
    i = Ix[0]
    if U_0 is None:
        # Запись граничных условий (x=0: i-1 -> i+1  u[i-1]=u[i+1]
        # где du/dn = 0, on x=L: i+1 -> i-1  u[i+1]=u[i-1])
        ip1 = i + 1
        im1 = ip1  # i-1 -> i+1
        u[i] = u_n[i] + dt * V(x[i]) + \
               0.5 * C2 * (0.5 * (q[i] + q[ip1]) * (u_n[ip1] - u_n[i]) - 0.5 * (q[i] + q[im1])
                           * (u_n[i] - u_n[im1])) + 0.5 * dt2 * f(x[i], t[0])
    else:
        u[i] = U_0(dt)
 
    i = Ix[-1]
    if U_L is None:
        im1 = i - 1
        ip1 = im1  # i+1 -> i-1
        u[i] = u_n[i] + dt * V(x[i]) + \
               0.5 * C2 * (0.5 * (q[i] + q[ip1]) * (u_n[ip1] - u_n[i]) - 0.5 * (q[i] + q[im1]) * (u_n[i] - u_n[im1])) + \
               0.5 * dt2 * f(x[i], t[0])
    else:
        u[i] = U_L(dt)
 
    if user_action is not None:
        user_func(u, x, t, 1)
 
    # Обновление данных и подготовка к новому шагу
    u_nm, u_n, u = u_n, u, u_nm
 
    # --- Симуляция (цикл прохода по времени) ---
    for n in It[1:-1]:
        # Обновление значений во внутренних узлах сетки
        for i in Ix[1:-1]:
            u[i] = - u_nm[i] + 2 * u_n[i] + \
                   C2 * (0.5 * (q[i] + q[i + 1]) * (u_n[i + 1] - u_n[i]) -
                         0.5 * (q[i] + q[i - 1]) * (u_n[i] - u_n[i - 1])) + dt2 * f(x[i], t[n])
 
        #  --- Запись граничных условий ---
        i = Ix[0]
        if U_0 is None:
            # Установка значений граничных условий
            # x=0: i-1 -> i+1  u[i-1]=u[i+1] где du/dn=0
            # x=L: i+1 -> i-1  u[i+1]=u[i-1] где du/dn=0
            ip1 = i + 1
            im1 = ip1
            u[i] = - u_nm[i] + 2 * u_n[i] + \
                   C2 * (0.5 * (q[i] + q[ip1]) * (u_n[ip1] - u_n[i]) - 0.5 * (q[i] + q[im1]) * (u_n[i] - u_n[im1])) + \
                   dt2 * f(x[i], t[n])
        else:
            u[i] = U_0(t[n + 1])
 
        i = Ix[-1]
        if U_L is None:
            im1 = i - 1
            ip1 = im1
            u[i] = - u_nm[i] + 2 * u_n[i] + \
                   C2 * (0.5 * (q[i] + q[ip1]) * (u_n[ip1] - u_n[i]) - 0.5 * (q[i] + q[im1]) * (u_n[i] - u_n[im1])) + \
                   dt2 * f(x[i], t[n])
        else:
            u[i] = U_L(t[n + 1])
 
        if user_action is not None:
            if user_action(u, x, t, n + 1):
                break
 
        u_nm, u_n, u = u_n, u, u_nm
 
    return u, x, t
 
 
# Функция формирования I(x)
def I(x):
    """
    Задается выражение для начального условия
    :param x:
    :return:
    """
 
    return 0.2 * (1 - x) * sin(math.pi * x)
 
 
# Функция формирования  V(x)
def V(x):
    """
    Задачется выражение для начальной производной
    :param x:
    :return:
    """
 
    return 0  # задается необходимое выражение
 
 
# Функция формирования  f(x,t)
def f(x, t):
    """
    Задается функция f(x,t), если данное уравнение не является однородным
    :param x:
    :param t:
    :return:
    """
 
    return 0  # задается необходимое выражение
 
 
# Функция для задания граничных условий
def U_0( ):
    """
    Задается граничное значение или функция (условие слева)
    :return:
    """
 
    return 0  # задается необходимое выражение
 
 
# Функция для задания граничных условий
def U_L( ):
    """
    Задается граничное значение или функция (условие справа)
    :return:
    """
 
    return 0  # задается необходимое выражение
 
 
# Функция запуска процесса симуляции
def simulate(
        I, V, f, c, L, U_0, U_L, dt, C, T,  # Параметры модели (уравнения)
        umin, umax,  # Интервал амплитуды колебаний
        animate = True,  # Выбор опции визуализации (с анимацией или нет)
        solver_func = sol,  # Вызов решателя
        mode = 'plotter',  # Выбор режима работы: вывод графика или сохранение данных в файл
):
    """
    Функция запуска процесса симуляции и вывода графика численного решения
    волнового уравнения, а также записи результатов в файл.
    :param I:
    :param V:
    :param f:
    :param c:
    :param L:
    :param dt:
    :param C:
    :param T:
    :param umin:
    :param umax:
    :param U_0:
    :param U_L:
    :param animate:
    :param solver_func:
    :param mode:
    :return:
    """
 
    class PlotMatplotlib:
        def __call__(self, u, x, t, n):
            """
            Функция типа user_func для визуализации
            :param u:
            :param x:
            :param t:
            :param n:
            :return:
            """
            if n == 0:
                plt.ion()
                self.lines = plt.plot(x, u, 'r-')
                plt.xlabel('x')
                plt.ylabel('u')
                plt.axis([0, L, umin, umax])
                plt.legend(['t=%f' % t[n]], loc = 'lower left')
            else:
                self.lines[0].set_ydata(u)
                plt.legend(['t=%f' % t[n]], loc = 'lower left')
                plt.draw()
 
            # Начало формирования кадров (первый кадр)
            time.sleep(2) if t[n] == 0 else time.sleep(0.2)
            plt.savefig('tmp_%04d.png' % n)  # для создания анимации
 
    if mode == 'plotter':
        plot_u = PlotMatplotlib()
    elif mode == 'logger':
        ff = file_save()
 
    # Очистка фреймов для визуализации
    for filename in glob.glob('tmp_*.png'):
        os.remove(filename)
 
    # Вызов функции решателя и начало симуляции
    user_func = plot_u if animate else None
    u, x, t = solver_function(
            I, V, f, c, L, dt, U_0, U_L, C, T, user_func)
 
    # Создание видеофайла анимации
    fps = 5
 
    # перечень поддерживаемых видеокодеков
    codec2ext = {'flv': 'flv', 'libx264': 'mp4', 'libvpx': 'webm', 'libtheora': 'ogg'}
    filespec = 'tmp_%04d.png'
    movie_program = 'ffmpeg'  # или 'avconv'
    for codec in codec2ext:
        ext = codec2ext[codec]
        cmd = '{movie_program} -r {fps:d} -i {filespec} ' \
              '-vcodec {codec} movie.{ext}'.format(**vars())
        os.system(cmd)
 
    # Запись данных симуляции в файл
    def file_save( ):
        """
        Сохранение данных расечтов  в файл
        :return:
        """
 
        sys.stdout = open('data.txt', 'w')
        u, x, t = solver_func(
                I, V, f, c, L, U_0, U_L, dt, C, T)
        print(x, u, t)
        sys.stdout.close()
        return print('writing data is success')
 
    return 0
 
 
def task( ):
    '''
    Постановка задачи
    :return:
    '''
 
    L = 1
    a = 1
    C = 0.75
    T = 1
    U_0(), U_L(), V, f, I
    umax = 2
    umin = -umax
 
    simulate(I, V, f, a, L, dt, U_0, U_L, C, T, umax, umin, user_func,
             animate = True, mode = 'plotter', solver_func = 'sol')
 
 
if __name__ == '__main__':
    task()

Выдаваемая ошибка:

Traceback (most recent call last):
File "C:\\LR2-rep\wave_eq_1d.py", line 350, in <module>
task()
File "C:\\LR2-rep\wave_eq_1d.py", line 346, in task
simulate(I, V, f, a, L, dt, U_0, U_L, C, T, umax, umin, animate = True, mode = 'plotter', solver_func = 'sol')
File "C:\\LR2-rep\wave_eq_1d.py", line 299, in simulate
solver_func(I, V, f, a, L, dt, U_0, U_L, C, T, user_func)
TypeError: 'str' object is not callable

Заранее благодарю!

@u235 · 26.10.2021, 06:08

solver_func у тебя просто строка, а ты, видимо, пытаешься ее вызывать как функцию.

@dux99 · 26.10.2021, 09:34 **[ТС]**

u235,

Да, просто ночью уже совем не понимаешь, что происходит, прочитал Ваш ответ, надо было кавычки убрать где вызов sol был, но теперь я не понимаю другого: я обозначил переменную dx вроде бы правильно, и не понимаю, как он расценил параметр C на вход функции еще одну функцию, а не например число... Я же ведь присвоил С воплне конкретное определенное значение...

File "C:\\LR2-rep\wave_eq_1d.py", line 61, in sol
dx = dt * a / float(C)
TypeError: float() argument must be a string or a number, not 'function'

Добавлено через 25 минут
Но теперь я не понимаю другого: теперь ошибка получилась еще более прозаичная, не понимаю ее...?

File "C:\\LR2-rep\wave_eq_1d.py", line 126, in sol
u[i] = U_L(dt)
TypeError: 'float' object is not callable

@dux99 · 27.10.2021, 18:45 **[ТС]**

Заметил еще несколько ошибок, в том числе неправильный порядок аргументов, вот правленный код, однако теперь код, правда появилась еще одна ошибка...

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
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
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
# Author: Alm99-collab
# Date: 25/10/2021 2:00
# Description: MSUT STANKIN Μyagkov Alexandr IDM - 21 - 03
 
 
"""
Программа для численного решения волнового уравнения в одномерном случае,
при разлинчых начальных и граничных условий, с помощью явной разностной
схемы типа "крест".
 
Общий вид волнового уравнения:
u_tt = a**2*u_xx + f(x,t) (0,L) где u=0 на диапазоне x=0,L, for t in (0,T].
Начальные условия в общем случае: u=I(x), u_t=V(x).
В случае неоднородного уравнения задается функция f(x,t).
 
Синтакис основной функции решателя:
u, x, t = sol(I, V, f, a, L, dt, C, T, user_func) где:
I = I(x) - функция.
V = V(x) - функция.
f = f(x,t) - функция.
U_0, U_L, - условия на границе.
C - число Куранта (=c*dt/dx), зависящее от шага dx. Является криетрием
стабильности численных расчето, если соблюдено условие (<=1)
dt - шаг по времени.
dx - шаг по координате.
T - конечное время симуляции волнового процесса.
user_func - функция (u, x, t, n) для вызова пользовательских сценариев,
таких как анимация или вывод графика, запси данных в текстовый файл,
расчет ошибки (в случае если известно точное решение) и.т.д.
"""
 
import numpy as np
import math
import matplotlib.pyplot as plt
import os
import time
import glob
 
 
def sol(I, V, f, a, L, C, T, U_0, U_L, dt, user_func = None):
    """
    Функция решатель волнового 1-D волнового уравнения
    u_tt = a**2*u_xx + f(x,t) (0,L) где u=0 на диапазоне
    x=0,L, for t in (0,T].
    :param I:
    :param V:
    :param f:
    :param a:
    :param L:
    :param C:
    :param T:
    :param U_0:
    :param U_L:
    :param dt:
    :param user_func:
    :return:
    """
 
    nt = int(round(T / dt))
    t = np.linspace(0, nt * dt, nt + 1)  # Массив точек для сетки времени
    dx = dt * a / float(C)
    nx = int(round(L / dx))
    x = np.linspace(0, L, nx + 1)  # Массив точек для сетки координат
    q = a ** 2
    C2 = (dt / dx) ** 2
    dt2 = dt * dt
 
    # --- Условие для проверки в качестве входного аргумента однородной функии ---
    if f is None or f == 0:
        f = lambda x, t: 0  # без использования векторизации используем лямбда- функцию
 
    # --- Условие для проверки в качестве входного аргумента начального условия dU(x,0)/dt ---
    if V is None or V == 0:
        V = lambda x: 0
 
    # Граничные условия для заданного волнового уравнения
    if U_0 is not None:
        if isinstance(U_0, (float, int)) and U_0 == 0:
            U_0 = lambda t: 0
    if U_L is not None:
        if isinstance(U_L, (float, int)) and U_L == 0:
            U_L = lambda t: 0
 
    # ---  Выделяем память под значения решений  ---
    u = np.zeros(nx + 1)  # Массив решений в узлах сетки на  временном шаге u(i,n+1)
    u_n = np.zeros(nx + 1)  # Массив решений в узлах сетки на  временном шаге u(i,n)
    u_nm = np.zeros(nx + 1)  # Массив решений в узлах сетки на  временном шаге u(i,n-1)
 
    # --- Проверка индексов для соблюдения размерностей массивов ---
    Ix = range(0, nx + 1)
    It = range(0, nt + 1)
 
    # --- Запись начальных условий ---
    for i in range(0, nx + 1):
        u_n[i] = I(x[i])
 
    if user_func is not None:
        user_func(u_n, x, t, 0)
 
    # --- Разностная формулма явной схемы "типа крест" на первом шаге ---
    for i in Ix[1:-1]:
        u[i] = u_n[i] + dt * V(x[i]) + 0.5 * C2 * (0.5 * (q[i] + q[i + 1]) * (u_n[i + 1] - u_n[i]) -
               0.5 * (q[i] + q[i - 1]) * (u_n[i] - u_n[i - 1])) + 0.5 * dt2 * f(x[i], t[0])
 
    i = Ix[0]
    if U_0 is None:
        # Запись граничных условий (x=0: i-1 -> i+1  u[i-1]=u[i+1]
        # где du/dn = 0, on x=L: i+1 -> i-1  u[i+1]=u[i-1])
        ip1 = i + 1
        im1 = ip1  # i-1 -> i+1
        u[i] = u_n[i] + dt * V(x[i]) + \
               0.5 * C2 * (0.5 * (q[i] + q[ip1]) * (u_n[ip1] - u_n[i]) - 0.5 * (q[i] + q[im1])
                           * (u_n[i] - u_n[im1])) + 0.5 * dt2 * f(x[i], t[0])
    else:
        u[i] = U_0(dt)
 
    i = Ix[-1]
    if U_L is None:
        im1 = i - 1
        ip1 = im1  # i+1 -> i-1
        u[i] = u_n[i] + dt * V(x[i]) + \
               0.5 * C2 * (0.5 * (q[i] + q[ip1]) * (u_n[ip1] - u_n[i]) - 0.5 * (q[i] + q[im1]) * (u_n[i] - u_n[im1])) + \
               0.5 * dt2 * f(x[i], t[0])
    else:
        u[i] = U_L(dt)
 
    if user_func is not None:
        user_func(u, x, t, 1)
 
    # Обновление данных и подготовка к новому шагу
    u_nm, u_n, u = u_n, u, u_nm
 
    # --- Симуляция (цикл прохода по времени) ---
    for n in It[1:-1]:
        # Обновление значений во внутренних узлах сетки
        for i in Ix[1:-1]:
            u[i] = - u_nm[i] + 2 * u_n[i] + \
                   C2 * (0.5 * (q[i] + q[i + 1]) * (u_n[i + 1] - u_n[i]) -
                         0.5 * (q[i] + q[i - 1]) * (u_n[i] - u_n[i - 1])) + dt2 * f(x[i], t[n])
 
        #  --- Запись граничных условий ---
        i = Ix[0]
        if U_0 is None:
            # Установка значений граничных условий
            # x=0: i-1 -> i+1  u[i-1]=u[i+1] где du/dn=0
            # x=L: i+1 -> i-1  u[i+1]=u[i-1] где du/dn=0
            ip1 = i + 1
            im1 = ip1
            u[i] = - u_nm[i] + 2 * u_n[i] + \
                   C2 * (0.5 * (q[i] + q[ip1]) * (u_n[ip1] - u_n[i]) - 0.5 * (q[i] + q[im1]) * (u_n[i] - u_n[im1])) + \
                   dt2 * f(x[i], t[n])
        else:
            u[i] = U_0(t[n + 1])
 
        i = Ix[-1]
        if U_L is None:
            im1 = i - 1
            ip1 = im1
            u[i] = - u_nm[i] + 2 * u_n[i] + \
                   C2 * (0.5 * (q[i] + q[ip1]) * (u_n[ip1] - u_n[i]) - 0.5 * (q[i] + q[im1]) * (u_n[i] - u_n[im1])) + \
                   dt2 * f(x[i], t[n])
        else:
            u[i] = U_L(t[n + 1])
 
        if user_func is not None:
            if user_func(u, x, t, n + 1):
                break
 
        u_nm, u_n, u = u_n, u, u_nm
 
    return u, x, t
 
 
# Функция формирования I(x)
def I(x):
    """
    Задается выражение для начального условия
    :param x:
    :return:
    """
 
    return 0.2 * (1 - x) * math.sin(math.pi * x)
 
 
# Функция формирования  V(x)
def V(x):
    """
    Задачется выражение для начальной производной
    :param x:
    :return:
    """
 
    return 0  # задается необходимое выражение
 
 
# Функция формирования  f(x,t)
def f(x, t):
    """
    Задается функция f(x,t), если данное уравнение не является однородным
    :param x:
    :param t:
    :return:
    """
 
    return 0  # задается необходимое выражение
 
 
# Функция для задания граничных условий
def U_0( ):
    """
    Задается граничное значение или функция (условие слева)
    :return:
    """
 
    return 0  # задается необходимое выражение
 
 
# Функция для задания граничных условий
def U_L( ):
    """
    Задается граничное значение или функция (условие справа)
    :return:
    """
 
    return 0  # задается необходимое выражение
 
 
# Функция запуска процесса симуляции
def simulate(
        I, V, f, a, L, C, T, U_0, U_L, dt,  # Параметры модели (уравнения)
        umin, umax,  # Интервал амплитуды колебаний
        animate = True,  # Выбор опции визуализации (с анимацией или нет)
        solver_func = sol,  # Вызов решателя
        mode = 'plotter',  # Выбор режима работы: вывод графика или сохранение данных в файл
):
    """
    Функция запуска процесса симуляции и вывода графика численного решения
    волнового уравнения, а также записи результатов в файл.
    :param I:
    :param V:
    :param f:
    :param a:
    :param L:
    :param dt:
    :param C:
    :param T:
    :param umin:
    :param umax:
    :param U_0:
    :param U_L:
    :param animate:
    :param solver_func:
    :param mode:
    :return:
    """
 
    class PlotMatplotlib:
        def __call__(self, u, x, t, n):
            """
            Функция типа user_func для визуализации
            :param u:
            :param x:
            :param t:
            :param n:
            :return:
            """
            if n == 0:
                plt.ion()
                self.lines = plt.plot(x, u, 'r-')
                plt.xlabel('x')
                plt.ylabel('u')
                plt.axis([0, L, umin, umax])
                plt.legend(['t=%f' % t[n]], loc = 'lower left')
            else:
                self.lines[0].set_ydata(u)
                plt.legend(['t=%f' % t[n]], loc = 'lower left')
                plt.draw()
 
            # Начало формирования кадров (первый кадр)
            time.sleep(2) if t[n] == 0 else time.sleep(0.2)
            plt.savefig('tmp_%04d.png' % n)  # для создания анимации
 
    if mode == 'plotter':
        plot_u = PlotMatplotlib()
    elif mode == 'logger':
        ff = file_save()
 
    # Очистка фреймов для визуализации
    for filename in glob.glob('tmp_*.png'):
        os.remove(filename)
 
    # Вызов функции решателя и начало симуляции
    user_func = plot_u if animate else None
    solver_func(I, V, f, a, L, C, T, U_0, U_L, dt, user_func)
 
    # Создание видеофайла анимации
    fps = 5
 
    # перечень поддерживаемых видеокодеков
    codec2ext = {'flv': 'flv', 'libx264': 'mp4', 'libvpx': 'webm', 'libtheora': 'ogg'}
    filespec = 'tmp_%04d.png'
    movie_program = 'ffmpeg'  # или 'avconv'
    for codec in codec2ext:
        ext = codec2ext[codec]
        cmd = '{movie_program} -r {fps:d} -i {filespec} ' \
              '-vcodec {codec} movie.{ext}'.format(**vars())
        os.system(cmd)
 
    # Запись данных симуляции в файл
    def file_save( ):
        """
        Сохранение данных расечтов  в файл
        :return:
        """
 
        sys.stdout = open('data.txt', 'w')
        u, x, t = solver_func(
                I, V, f, a, L, C, T, U_0, U_L, dt)
        print(x, u, t)
        sys.stdout.close()
        return print('writing data is success')
 
    return 0
 
 
def task( ):
    '''
    Постановка задачи
    :return:
    '''
 
    I
    L = 1
    a = 1
    C = 0.85
    T = 1
    dt = 0.05
    U_0, U_L, V, f = 0, 0, 0, 0
    umax = 2
    umin = -umax
 
    simulate(I, V, f, a, L, C, T, U_0, U_L, dt, umax, umin, animate = True, solver_func = sol, mode = 'plotter')
 
 
if __name__ == '__main__':
    task()

Добавлено через 55 секунд
u235,

После последней правки ошибка такая:

Code
1
2
3
  File "C:\\LR2-rep\wave_eq_1d.py", line 102, in sol
    u[i] = u_n[i] + dt * V(x[i]) + 0.5 * C2 * (0.5 * (q[i] + q[i + 1]) * (u_n[i + 1] - u_n[i]) -
TypeError: 'int' object is not subscriptable

Добавлено через 27 минут
u235,

Ошибка на строке 338, при переопределении U_0 и U_L как чисел, хотя они являются функциями. Должны ими быть.

В данной постановке, граничные условия должны быть равны нулю, но в задании также сказано провести расчеты при произвольных параметрах также, т.е. исследовать влияние граничных и начальных условий, в том числе если они будут заданы в виде функций... Или мне тогда просто лучше в самой функции task посредством лямбда-выражений определять их и передавать?

Как бы эту программу я хочу использовать еще и для решения другой задачи, где приводится тоже волнвое уравнение второго порядка, однако, там добавлена еще функция и граничные условия заданы на одном конце в виде ступенчатой функции.

Добавлено через 4 часа 54 минуты
Я переписал функцию task():

Python
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
def task( ):
    '''
    Постановка задачи
    :return:
    '''
 
    I
    L = 1
    a = 1
    C = 0.85
    T = 1
    dt = 0.05
    U_0, U_L, V, f
    umax = 2
    umin = -umax
 
    simulate(I, V, f, a, L, C, T, U_0, U_L, dt, umax, umin, animate = True, solver_func = sol, mode = 'plotter',)
 
 
if __name__ == '__main__':
    task()

И получил ошибку следующего типа:

Code
1
2
3
File "C:\\LR2-rep\wave_eq_1d.py", line 102, in sol
    u[i] = u_n[i] + dt * V(x[i]) + 0.5 * C2 * (0.5 * (q[i] + q[i + 1]) * (u_n[i + 1] - u_n[i]) -
TypeError: 'int' object is not subscriptable

Ошибка в принципе мне сама понятна, что я пытаюсь рассматривать объекти типа int как например объект типов list,str, map... Меня просто интересует, ведь в принципе я просто писал шаблон для разностной схемы для элментов , которые индексируеются, как например в matlab, и вроде бы правильно делаю, но получается, что мне придется задавать u_nm, u как элменты например list? Почему-то считает, что это объект целочисленный, и не поддается подписи... Когда мы пытаемся объединить строковые и целочисленные значения, это сообщение говорит нам, что мы рассматриваем целое число как подписываемый объект. Целое число не является подписываемым объектом, но как тогда еще исправлять подобное?

@dux99 · 29.10.2021, 21:11 **[ТС]**

u235,

В общем попытался переписать код иначе, единственный вопрос у меня возникает в том,что я не уверен, что решение получено точное, и еще как бы вывести в формате столбцов ndarray? К сожалению не получается из отдельных кадров создать анимацию, к тому же не на всех кадрах прорисовывается сетка... Не знаю как это исправить.... Может есть какие предложения?

Прикладываю код программы, потому что хотелось бы и иметь видеоанимацию, так и по верменные фреймы, для отслеживания решения...

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
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
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
# Author: Alm99-collab
# Date: 25/10/2021 2:00
# Description: MSUT STANKIN Μyagkov Alexandr IDM - 21 - 03
 
 
"""
Программа для численного решения волнового уравнения в одномерном случае,
при разлинчых начальных и граничных условий, с помощью явной разностной
схемы типа "крест".
 
Общий вид волнового уравнения:
u_tt = a**2*u_xx + f(x,t) (0,L) где u=0 на диапазоне x=0,L, for t in (0,T].
Начальные условия в общем случае: u=I(x), u_t=V(x).
В случае неоднородного уравнения задается функция f(x,t).
 
Синтакис основной функции решателя:
u, x, t = solver(I, V, f, a, L, dt, C, T, user_func) где:
I = I(x) - функция.
V = V(x) - функция.
f = f(x,t) - функция.
U_0, U_L, - условия на границе.
C - число Куранта (=a*dt/dx), зависящее от шага dx. Является криетрием
стабильности численных расчето, если соблюдено условие (<=1)
dt - шаг по времени.
dx - шаг по координате.
T - конечное время симуляции волнового процесса.
user_func - функция (u, x, t, n) для вызова пользовательских сценариев,
таких как анимация или вывод графика, запси данных в текстовый файл,
расчет ошибки (в случае если известно точное решение) и.т.д.
"""
 
import numpy as np
import math
import matplotlib.pyplot as plt
import os
import time
import glob
 
 
# функция решатель для данного дифференциального уравнения
def solver(I, V, f, a, U_0, U_L, L, dt, C, T,
           user_func = None):
    """
    Функция решения волнового уравнения u_tt = a**2*u_xx + f(x,t) (0,L),
    где u=0 на диапазоне x=0,L, for t in (0,T].
    Начальные условия в общем случае: u=I(x), u_t=V(x).
    В случае неоднородного уравнения задается функция f(x,t).
    ------------------------------------------------------------------------
    :param I:
    :param V:
    :param f:
    :param a:
    :param L:
    :param C:
    :param T:
    :param U_0:
    :param U_L:
    :param dt:
    :param user_func:
    :return:
    """
 
    nt = int(round(T / dt))
    t = np.linspace(0, nt * dt, nt + 1)  # Узлы сетки по времени
    dx = dt * a / float(C)
    nx = int(round(L / dx))
    x = np.linspace(0, L, nx + 1)  # Узлы сетки по координате
    C2 = C ** 2
    dt2 = dt * dt
 
    # Проверка того, что  массивы являются элементами t,x
    dx = x[1] - x[0]
    dt = t[1] - t[0]
 
    # Выбор и инициализация дополнительных параметров f, I, V, U_0, U_L если равны нулю
    # или не передаются
    if f is None or f == 0:
        f = lambda x, t: 0
    if I is None or I == 0:
        I = lambda x: 0
    if V is None or V == 0:
        V = lambda x: 0
    if U_0 is not None:
        if isinstance(U_0, (float, int)) and U_0 == 0:
            U_0 = lambda t: 0
        # иначе: U_0(t) является функцией
    if U_L is not None:
        if isinstance(U_L, (float, int)) and U_L == 0:
            U_L = lambda t: 0
        # иначе: U_L(t) является функцией
 
    # ---  Выделяем память под значения решений  ---
    u = np.zeros(nx + 1)  # Массив решений в узлах сетки на  временном шаге u(i,n+1)
    u_n = np.zeros(nx + 1)  # Массив решений в узлах сетки на  временном шаге u(i,n)
    u_nm1 = np.zeros(nx + 1)  # Массив решений в узлах сетки на  временном шаге u(i,n-1)
 
    # --- Проверка индексов для соблюдения размерностей массивов ---
    Ix = range(0, nx + 1)
    It = range(0, nt + 1)
 
    # --- Запись начальных условий ---
    for i in Ix:
        u_n[i] = I(x[i])
 
    if user_func is not None:
        user_func(u_n, x, t, 0)
 
    # --- Разностная формулма явной схемы "типа крест" на первом шаге ---
    for i in Ix[1:-1]:
        u[i] = u_n[i] + dt * V(x[i]) + 0.5 * C2 * (u_n[i - 1] - 2 * u_n[i] + u_n[i + 1]) + 0.5 * dt2 * f(x[i], t[0])
 
    i = Ix[0]
    if U_0 is None:
        # Запись граничных условий (x=0: i-1 -> i+1  u[i-1]=u[i+1]
        # где du/dn = 0, on x=L: i+1 -> i-1  u[i+1]=u[i-1])
        ip1 = i + 1
        im1 = ip1  # i-1 -> i+1
        u[i] = u_n[i] + dt * V(x[i]) + 0.5 * C2 * (u_n[im1] - 2 * u_n[i] + u_n[ip1]) + 0.5 * dt2 * f(x[i], t[0])
 
    else:
        u[0] = U_0(dt)
 
    i = Ix[-1]
    if U_L is None:
        im1 = i - 1
        ip1 = im1  # i+1 -> i-1
        u[i] = u_n[i] + dt * V(x[i]) + 0.5 * C2 * (u_n[im1] - 2 * u_n[i] + u_n[ip1]) + 0.5 * dt2 * f(x[i], t[0])
    else:
        u[i] = U_L(dt)
 
    if user_func is not None:
        user_func(u_n, x, t, 1)
 
    # Обновление данных и подготовка к новому шагу
    u_nm1, u_n, u = u_n, u, u_nm1
 
    # --- Симуляция (цикл прохода по времени) ---
    for n in It[1:-1]:
        # Обновление значений во внутренних узлах сетки
        for i in Ix[1:-1]:
            u[i] = - u_nm1[i] + 2 * u_n[i] + C2 * (u_n[i - 1] - 2 * u_n[i] + u_n[i + 1]) + dt2 * f(x[i], t[n])
 
        #  --- Запись граничных условий ---
        i = Ix[0]
        if U_0 is None:
            # Установка значений граничных условий
            # x=0: i-1 -> i+1  u[i-1]=u[i+1] где du/dn=0
            # x=L: i+1 -> i-1  u[i+1]=u[i-1] где du/dn=0
            ip1 = i + 1
            im1 = ip1
            u[i] = - u_nm1[i] + 2 * u_n[i] + C2 * (u_n[im1] - 2 * u_n[i] + u_n[ip1]) + dt2 * f(x[i], t[n])
        else:
            u[0] = U_0(t[n + 1])
 
        i = Ix[-1]
        if U_L is None:
            im1 = i - 1
            ip1 = im1
            u[i] = - u_nm1[i] + 2 * u_n[i] + C2 * (u_n[im1] - 2 * u_n[i] + u_n[ip1]) + dt2 * f(x[i], t[n])
        else:
            u[i] = U_L(t[n + 1])
 
        if user_func is not None:
            if user_func(u, x, t, n + 1):
                break
 
        # Обновление данных и подготовка к новому шагу
        u_nm1, u_n, u = u_n, u, u_nm1
 
    # Присвоение значений требуемому узлу после прохода по сетке
    u = u_n
 
    return u, x, t
 
 
# функция симуляции и анимации, сохранения данных в файл
def simulate(I, V, f, a, U_0, U_L, L, dt, C, T, umin, umax, animate = True):
 
    """
    Запуск решателя и анимации,сохранения данных в файл
    ----------------------------------------------------------------------------
    :param I:
    :param V:
    :param f:
    :param a:
    :param L:
    :param dt:
    :param C:
    :param T:
    :param umin:
    :param umax:
    :param U_0:
    :param U_L:
    :param animate:
    :param mode:      # выбор режима работы (анимация и визуализация или запись данных в файл)
    :return:
    """
 
    if callable(U_0):
        bc_left = 'u(0,t)=U_0(t)'
    elif U_0 is None:
        bc_left = 'du(0,t)/dx=0'
    else:
        bc_left = 'u(0,t)=0'
    if callable(U_L):
        bc_right = 'u(L,t)=U_L(t)'
    elif U_L is None:
        bc_right = 'du(L,t)/dx=0'
    else:
        bc_right = 'u(L,t)=0'
 
    class PlotMatplotlib:
        def __call__(self, u, x, t, n):
 
            """ user_func для визуализации """
 
            if n == 0:
                plt.ion()
                self.lines = plt.plot(x, u, 'r-')
                plt.xlabel('x')
                plt.ylabel('u')
                plt.axis([0, L, umin, umax])
                plt.legend(['t=%f' % t[n]], loc = 'lower left')
            else:
                self.lines[0].set_ydata(u)
                plt.legend(['t=%f' % t[n]], loc = 'lower left')
                plt.draw()
                plt.grid()
            time.sleep(1) if t[n] == 0 else time.sleep(0.5)
            plt.savefig('tmp_%04d.png' % n)  # для создания анимации из кадров
 
    plot_u = PlotMatplotlib()
 
    # Очистка фреймов для создания анимированного изображения
    for filename in glob.glob('frame_*.png'):
        os.remove(filename)
 
    fps = 4  # frames per second
    codec2ext = dict(flv = 'flv', libx264 = 'mp4', libvpx = 'webm',
                     libtheora = 'ogg')  # video formats
    filespec = 'tmp_%04d.png'
    movie_program = 'ffmpeg'  # or 'avconv'
    for codec in codec2ext:
        ext = codec2ext[codec]
        cmd = '%(movie_program)s -r %(fps)d -i %(filespec)s ' \
              '-vcodec %(codec)s movie.%(ext)s' % vars()
        os.system(cmd)
 
    user_func = plot_u if animate else None
    u, x, t = solver(I, V, f, a, U_0, U_L, L, dt, C, T, user_func)
 
    return u, x, t
 
 
# функция постановки задачи
def problem():
    I = lambda x: 0.2 * (1 - x) * math.sin(math.pi * x)
    V = lambda x: 0
    f = lambda x, t: 0
    U_0 = lambda t: 0
    U_L = 0
    L = 1
    a = 1
    C = 0.75
    nx = 12
    dt = C * (L / nx) / a
    T = 1
    umin = -0.25
    umax = 0.25
 
    u, x, t = simulate(I, V, f, a, U_0, U_L, L, dt, C, T, umin, umax, animate=True)
 
    return u, x, t
 
 
if __name__ == '__main__':
    u, x, t = problem()
    print(u,x,t)

@dux99 · 31.10.2021, 15:49 **[ТС]**

В целом, удалось решить задачу, и создать анимацию!

Новые блоги и статьи Все статьи Все блоги /
PhpStorm 2025.3: WSL Terminal всегда стартует в ~ and_y87 14.12.2025 PhpStorm 2025. 3: WSL Terminal всегда стартует в ~ (home), игнорируя директорию проекта Симптом: После обновления до PhpStorm 2025. 3 встроенный терминал WSL открывается в домашней директории. . .	Как объединить две одинаковые БД Access с разными данными VikBal 11.12.2025 Помогите пожалуйста !! Как объединить 2 одинаковые БД Access с разными данными.	Новый ноутбук volvo 07.12.2025 Всем привет. По скидке в "черную пятницу" взял себе новый ноутбук Lenovo ThinkBook 16 G7 на Амазоне: Ryzen 5 7533HS 64 Gb DDR5 1Tb NVMe 16" Full HD Display Win11 Pro	Музыка, написанная Искусственным Интеллектом volvo 04.12.2025 Всем привет. Некоторое время назад меня заинтересовало, что уже умеет ИИ в плане написания музыки для песен, и, собственно, исполнения этих самых песен. Стихов у нас много, уже вышли 4 книги, еще 3. . .	От async/await к виртуальным потокам в Python IndentationError 23.11.2025 Армин Ронахер поставил под сомнение async/ await. Создатель Flask заявляет: цветные функции - провал, виртуальные потоки - решение. Не threading-динозавры, а новое поколение лёгких потоков. Откат?. . .
Поиск "дружественных имён" СОМ портов Argus19 22.11.2025 Поиск "дружественных имён" СОМ портов На странице: https:/ / norseev. ru/ 2018/ 01/ 04/ comportlist_windows/ нашёл схожую тему. Там приведён код на С++, который показывает только имена СОМ портов, типа,. . .	Сколько Государство потратило денег на меня, обеспечивая инсулином. Programma_Boinc 20.11.2025 Сколько Государство потратило денег на меня, обеспечивая инсулином. Вот решила сделать интересный приблизительный подсчет, сколько государство потратило на меня денег на покупку инсулинов. . . .	Ломающие изменения в C#.NStar Alpha Etyuhibosecyu 20.11.2025 Уже можно не только тестировать, но и пользоваться C#. NStar - писать оконные приложения, содержащие надписи, кнопки, текстовые поля и даже изображения, например, моя игра "Три в ряд" написана на этом. . .	Мысли в слух kumehtar 18.11.2025 Кстати, совсем недавно имел разговор на тему медитаций с людьми. И обнаружил, что они вообще не понимают что такое медитация и зачем она нужна. Самые базовые вещи. Для них это - когда просто люди. . .	Создание Single Page Application на фреймах krapotkin 16.11.2025 Статья исключительно для начинающих. Подходы оригинальностью не блещут. В век Веб все очень привыкли к дизайну Single-Page-Application . Быстренько разберем подход "на фреймах". Мы делаем одну. . .

@u235 5514 / 2867 / 571 Регистрация: 07.11.2019 Сообщений: 4,751
	26.10.2021, 06:08
	solver_func у тебя просто строка, а ты, видимо, пытаешься ее вызывать как функцию. 1

@dux99 8 / 6 / 2 Регистрация: 03.05.2016 Сообщений: 242
	26.10.2021, 09:34 [ТС]
	u235, Да, просто ночью уже совем не понимаешь, что происходит, прочитал Ваш ответ, надо было кавычки убрать где вызов sol был, но теперь я не понимаю другого: я обозначил переменную dx вроде бы правильно, и не понимаю, как он расценил параметр C на вход функции еще одну функцию, а не например число... Я же ведь присвоил С воплне конкретное определенное значение... File "C:\\LR2-rep\wave_eq_1d.py", line 61, in sol dx = dt * a / float(C) TypeError: float() argument must be a string or a number, not 'function' Добавлено через 25 минут Но теперь я не понимаю другого: теперь ошибка получилась еще более прозаичная, не понимаю ее...? File "C:\\LR2-rep\wave_eq_1d.py", line 126, in sol u[i] = U_L(dt) TypeError: 'float' object is not callable 0

@dux99 8 / 6 / 2 Регистрация: 03.05.2016 Сообщений: 242
	31.10.2021, 15:49 [ТС]
	В целом, удалось решить задачу, и создать анимацию! 1