# 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()
