Попытка численно решить уравнение теплопроводности методом конечных разностей. Решение не удовлетворяет начальному условию. Есть идеи, почему?
Я подозреваю, что сюжет относится к самому начальному условию, а не к решению. Формула дискретизации верна.
Код
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import cm
from mpl_toolkits.mplot3d import Axes3D
M = 50 # number of grid points for space interval
N = 200 # '' '' '' '' '' time ''
x0 = 0
xL = 5 # grid differences
dx = (xL - x0) / (M - 1) # space step
t0 = 0
tF = 0.1 # grid differences
dt = (tF - t0) / (N - 1) # time step
D = 0.1 # thermal diffusivity
a = D * dt / dx**2
# Create grid
t = np.linspace(t0, tF, N)
x = np.linspace(x0, xL, M)
# Initial matrix solution
U = np.zeros((M, N))
# Initial condition
U[:, 0] = np.sin(np.pi * x)
# Boundary conditions
U[0, :] = 0
U[-1, :] = 0 # Dirichlet boundary conditions
# Discretised derivative formula
for k in range(0, N-1):
for i in range(1, M-1):
U[i, k+1] = a * U[i-1, k] + (1 - 2 * a) * U[i, k] + a * U[i + 1, k] # This is the formula itself!
#U[i - 1, k] = U[i - 3, k] # Von Neumann boundary condition
T, X = np.meshgrid(t, x)
fig = plt.figure()
ax = fig.gca(projection='3d')
soln = ax.plot_surface(X, T, U, cmap=cm.coolwarm,
linewidth=0, antialiased=False)
fig.colorbar(soln)
ax.set_xlabel('Time')
ax.set_ylabel('Space')
ax.set_zlabel('U(x, t) - Concentration')
plt.tight_layout()
plt.show()