Базовое одномерное уравнение адвекции в Python с использованием метода конечных разностей ⇐ Python

[Гость]

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

Код: Выделить всё

import numpy as np
import matplotlib.pyplot as plt

#Mesh
Nx = 211 #Number of x grid points
Nt = 1000 #Number of t grid points

dx = 1/(Nx - 1) #Dividing dx into equal dx intervales
dt = 1e-7 #Time step

#Boundary conditions
x= np.linspace(0.1, 1.0, num=Nx)
t= np.linspace(0.0, Nt*dt, num=Nt)

#Initializing array to zero
u=np.zeros((Nx, Nt)) #Initializing array to zero

#Velocity
c = 1

#Gaussian initial condition
mu = 0.5
sigma = 0.1

for i in range(1,Nx-1):
u[i,0] = np.exp(-(x[i]-mu)**2/(2*sigma**2)) / np.sqrt(2*np.pi*sigma**2)

#Looping over all grid points excluding those at the boundary conditions
for j in range(0,Nt-1):
for i in range(1,Nx-1):

u[i+1,j] = u[i,j] + c*(((1/x[i]) + (u[i+1,j] - (u[i,j])/dx)))*dt

plt.plot(x, u[:, 0])
plt.plot(x, u[:, 200])
plt.plot(x, u[:, 300])
plt.plot(x, u[:, 400])

Я ожидал, что график сдвинется влево, поскольку скорость c>0, и хотел увидеть график в 4 разных момента времени, начиная с начальной начальной точки по Гауссу. Но то, что я получаю, не имеет никакого смысла. Я также хочу отметить, что когда я комментирую цикл for и выполняю plt.plot(x,u[:,0] , я вижу начальную кривую Гаусса, как и ожидалось, когда я набираю ее снова после раскомментирования цикла for, я получаю график показано ниже, я был бы очень признателен за любую помощь, которую кто-либо может мне оказать в этом вопросе.
Выходной график

Подробнее здесь: https://stackoverflow.com/questions/784 ... nce-method