Решение PDE не соответствует начальному условию - PullRequest
0 голосов
/ 13 июля 2020

Попытка численно решить уравнение теплопроводности методом конечных разностей. Решение не удовлетворяет начальному условию. Есть идеи, почему?

Я подозреваю, что сюжет относится к самому начальному условию, а не к решению. Формула дискретизации верна.

Код

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