Дискуссионный математический форумМатематический форум
Математический форум Math Help Planet

Обсуждение и решение задач по математике, физике, химии, экономике

Теоретический раздел
Часовой пояс: UTC + 3 часа [ Летнее время ]
новый онлайн-сервис
число, сумма и дата прописью

Часовой пояс: UTC + 3 часа [ Летнее время ]




Начать новую тему Ответить на тему  [ Сообщений: 2 ] 
Автор Сообщение
 Заголовок сообщения: Уравнение переноса, явная схема
СообщениеДобавлено: 29 окт 2020, 22:02 
Не в сети
Одарённый
Зарегистрирован:
20 апр 2019, 13:04
Сообщений: 128
Cпасибо сказано: 19
Спасибо получено:
3 раз в 3 сообщениях
Очков репутации: 2

Добавить очки репутацииУменьшить очки репутации
Здравствуйте, исходная задача:

[math]\begin{array}{c}
\frac{\partial u}{\partial t}+a \frac{\partial u}{\partial x}=f(x, t), \quad a>0 \\
0 \leqslant x \leqslant 1,0 \leqslant t \leqslant 1 \\
u(x, 0)=\varphi(x), \quad x \in[0,1] \\
u(0, t)=g_{1}(t), \quad t \in[0,1]
\end{array}




\begin{array}{c}
\frac{U_{j}^{n+1}-U_{j}^{n}}{\tau}+a \frac{U_{3}^{n}-U_{j-1}^{n}}{h}=f_{j}^{n}, \\
j=\overline{1, M} ; \quad U_{0}^{n+1}=g_{1}^{n+1}
\end{array}



\begin{aligned}
a &=0.23 \\

u(x, t) &=x^{3}-\sin (2 \pi t x)*0.5+2.5 t x
\end{aligned}[/math]


Вот код на Python



Код:
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation

c = 0.23

xmin = 0
xmax = 1

n = 10

X, dx = np.linspace(xmin,xmax,n,retstep=True)
T, dt = np.linspace(xmin,xmax,n,retstep=True)



def u(x,t):
  return x*x*x - np.sin(2*np.pi*t*x)/2 + 2.5*t*x

def initial_u(x):
    return x*x*x

def f(x,t,a):
  return 2.5*x - np.pi*np.cos(2*np.pi*t*x) + a*(2*np.pi*np.pi*t*t*np.sin(2*np.pi*t*x) + 6*x)
 
U = []

# explicit euler solution
def u(x, t):
    if t == 0: # initial condition
        return initial_u(x)
    uvals = [] # u values for this time step
    for j in range(len(x)):
        if j == 0: # left boundary
            uvals.append(U[t-1][j] + c*dt/(2*dx)*(U[t-1][j+1]-U[t-1][n-1]))
        elif j == n-1: # right boundary
            uvals.append(U[t-1][j] + c*dt/(2*dx)*(U[t-1][0]-U[t-1][j-1]))
        else:
            uvals.append(U[t-1][j] + c*dt/(2*dx)*(U[t-1][j+1]-U[t-1][j-1]) + dt*f(X[j], T[t], c))
    return uvals

for t in range(10):
    U.append(u(X, t))
    plt.plot(X,U[-1])


Из разностной схемы выражаю:

[math]u_{j}^{n+1}=u_{j}^{n}+c \Delta t \frac{u_{j+1}^{n}-u_{j-1}^{n}}{2 \Delta x} + f(x,t)*\Delta t[/math]

Кажется все просто, но что-то не работает.. какая-то дичь получается. Что не так? Может кто-нибудь скинет рабочий код на любом языке? Я посмотрю, сравню.

Вернуться к началу
 Профиль  
Cпасибо сказано 
 Заголовок сообщения: Re: Уравнение переноса, явная схема
СообщениеДобавлено: 19 апр 2021, 17:00 
Не в сети
Гений
Зарегистрирован:
02 июн 2018, 08:50
Сообщений: 659
Cпасибо сказано: 21
Спасибо получено:
105 раз в 103 сообщениях
Очков репутации: 16

Добавить очки репутацииУменьшить очки репутации
Когда то писал решение уравнения теплопроводности с помощью МКР

[math]\frac{\partial U}{\partial t} = D \frac{\partial^2 U}{\partial x^2}[/math]

это не совсем то, но принцип похож

▼ код
module teplo
implicit none

public :: du_teplo
private

contains

subroutine du_teplo (d, nh, ns, t, u, x, root)
real :: q !Константа q
real :: d !Коэффициент теплопроводности материала
real :: h !Шаг по сетке в пространственной координате x
real :: s !Шаг по сетке во временной координате t
real :: t !Конечное время t
real :: temp !Временная переменная
real, dimension (3) :: u !Начальные значения температуры (1 - слева, 2 - справа, 3 - начальная по длине стержня)
real, dimension (2) :: x !Конечные значения пространственной координаты x (1-слева, 2 - справа)
real, dimension (nh) :: root !Массив значений температуры при конечном времени t
real, dimension (2,nh) :: z !Массив для вычислений
integer :: nh !Число шагов по сетке в пространственной координате
integer :: ns !Число шагов по сетке во временной координате t
integer :: i, j !Вспомогательные переменные

!Определяемся с исходными данными
z(1,:) = u(3)
z(:,1) = u(1)
z(:,nh) = u(2)
h = x(2) / float (nh)
s = t / float (ns)
temp = h*h / (2.0*d)

if (s >= temp) then !Обеспечим устойчивость
s = temp * 0.75
ns = nint(t/s)
end if

q = s / (h*h)

do j = 1, ns !Сделаем вычисления по сетке
do i = 2, nh-1
z(2,i) = q*d*z(1,i-1) + (1.0-2.0*d*q)*z(1,i) + q*d*z(1,i+1)
end do
z = cshift (z, 1)
end do

root = z(1,:)

end subroutine du_teplo
end module teplo

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

 Похожие темы   Автор   Ответы   Просмотры   Последнее сообщение 
Явная конечно-разностная схема

в форуме Численные методы

rokanten13

1

306

02 июн 2019, 14:06

Явная формула для предела f(x) в степени g(x)

в форуме Пределы числовых последовательностей и функций, Исследования функций

max_korostelev

1

140

27 сен 2021, 01:11

Программирование уравнения переноса

в форуме Информатика и Компьютерные науки

zolla

6

578

04 фев 2017, 17:57

Построение методом параллельного переноса

в форуме Геометрия

_Astarta_

5

711

28 апр 2018, 21:37

Решение уравнения переноса с помощью схемы Кранка-Николсона

в форуме Численные методы

Dellghin

2

954

03 мар 2017, 17:09

Как можно менять знаки переменных без переноса за знак равно

в форуме Алгебра

koderman

7

340

04 фев 2020, 08:17

Схема

в форуме Электричество и Магнетизм

DeWaldemar

0

361

07 июн 2015, 18:50

Схема

в форуме Дискретная математика, Теория множеств и Логика

drago123

4

292

28 ноя 2017, 18:00

Схема Бернулли?

в форуме Теория вероятностей

tanyhaftv

1

254

20 май 2018, 09:32

Схема Горнера

в форуме Информатика и Компьютерные науки

[Alexa]

0

762

03 июл 2021, 22:36


Часовой пояс: UTC + 3 часа [ Летнее время ]



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

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


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

Найти:
Перейти:  

Яндекс.Метрика

Copyright © 2010-2023 MathHelpPlanet.com. All rights reserved