在Python中使用fft求解热方程

Question

我正在尝试实现快速傅里叶变换来求解热方程。然而，该解决方案没有任何意义。

import numpy as np
from scipy import fft 
import matplotlib.pyplot as plt
from scipy.integrate import odeint
from mpl_toolkits.mplot3d import axes3d
import matplotlib.cm as cm
from scipy.integrate import solve_ivp

# Constants
c = 1  # thermal conductivity 
L = 4    # Length of the domain
N = 100  # Number of discretization points
dx = L / N
x = np.linspace(0, L, N, endpoint=False)

kappa = 2 * np.pi * fft.fftfreq(N, d=dx)


u0 = np.zeros_like(x)
start_idx = int(1 / dx)  # Index corresponding to x = 2
end_idx = int(3 / dx)   # Index corresponding to x = 10
u0[start_idx:end_idx] = 1  # Set the pulse between x = 2 and x = 10


def RHSheatspatial(t, u, kappa, c):
    uhat = fft.fft(u)
    dd_uhat = -np.power(kappa, 2) * uhat
    dd_u = fft.ifft(dd_uhat)
    du_dt = -c**2 * dd_u
    return np.real(du_dt)

t = np.linspace(0, 10, 1000) 

sol = solve_ivp(RHSheatspatial, (0, 10), u0, method='RK45', args=(kappa, c), t_eval=t,rtol=1.0e-6,atol=1.0e-8)
print(sol)

# Initial condition
plt.subplot(1, 2, 1)
plt.plot(x, u0, label='Initial Condition', color='blue')
plt.title('Initial Condition')
plt.xlabel('x')
plt.ylabel('T')
plt.grid(True)

# Final solution at t = 1
plt.subplot(1, 2, 2)
plt.plot(x, sol.y[:, -1], label='Final Solution at t=1', color='red')
plt.title('Solution at t=10')
plt.xlabel('x')
plt.ylabel('T')
plt.grid(True)

消息：所需的步长小于数字之间的间距。 成功：假 状态：-1 时间: [ 0.000e+00 1.001e-02 2.002e-02 3.003e-02 4.004e-02 5.005e-02 6.006e-02 7.007e-02 8.008e-02 9.009e-02 1.001e-01 1.101e-01] y: [[ 0.000e+00 -3.478e+22 ... -2.198e+254 -1.248e+280] [ 0.000e+00 1.042e+23 ... 6.584e+254 3.740e+280] ... [ 0.000e+00 1.042e+23 ... 6.584e+254 3.740e+280] [ 0.000e+00 -3.478e+22 ... -2.198e+254 -1.248e+280]] 溶胶：无 t_events：无 y_事件：无 NFEV：16994 新杰夫: 0 自然卢：0

Answer 1

您的

RHSheatspatial

中是否有过多的减号？当我删除其中一个时，我得到以下解决方案：

虽然我不确定这是否是你的热方程的正确解，但我相信它看起来更合理。