L2-норма может быть полезной метрикой для проверки сходимости, но не идеальна при отладке, поскольку не объясняет, в чем проблема.Хотя ваше решение должно быть безусловно стабильным, отсталый Эйлер не обязательно сходится к правильному ответу.Так же, как форвард Эйлер, как известно, нестабилен (антидиссипативен), отсталый Эйлер, как известно, диссипативен.Постановка ваших решений подтверждает это.Численные решения сходятся к нулю.Для приближения следующего порядка Крэнк-Николсон является разумным кандидатом.Приведенный ниже код содержит более общий тета-метод, так что вы можете настроить неявное решение.тета = 0,5 дает CN, тета = 1 дает BE, а тета = 0 дает FE.Несколько других вещей, которые я подправил:
- Я выбрал более подходящий временной шаг dt = (dx ** 2) / 2 вместо dt = dx.Это последнее не сходится к правильному решению с использованием CN.
- Это незначительное замечание, но поскольку t_final не гарантирует кратность dt, вы не сравнивали решения на одном шаге по времени.
- Что касается вашего комментария о том, что он медленный: по мере увеличения пространственного разрешения ваше разрешение по времени также должно увеличиваться.Даже в вашем случае с dt = dx вы должны выполнить (1024 x 1024) * 1024 матричное умножение 1024 раза.Я не нашел, что это займет особенно много времени на моей машине.Я удалил некоторую ненужную конкатенацию, чтобы немного ускорить ее, но изменение шага по времени на dt = (dx ** 2) / 2, к сожалению, действительно затормозит.Вы могли бы попытаться скомпилировать с Numba, если вы обеспокоены скоростью.
Все это говорит о том, что я не нашел огромного успеха с последовательностью CN.Я должен был установить N = 2 ^ 6, чтобы получить что-нибудь в t_final = 1.Увеличение t_final делает это хуже, уменьшение t_final делает его лучше.В зависимости от ваших потребностей, вы можете изучить реализацию TR-BDF2 или других линейных многошаговых методов, чтобы улучшить это.
Код с графиком приведен ниже:
import numpy as np
import matplotlib.pyplot as plt
tau = 2 * np.pi
tau2 = tau * tau
i = complex(0,1)
def solution_f(t, x):
return 0.5 * (np.exp(-tau * i * x) * np.exp((2 - tau2) * i * t) + np.exp(tau * i * x) * np.exp((tau2 + 4) * i * t))
def solution_g(t, x):
return 0.5 * (np.exp(-tau * i * x) * np.exp((2 - tau2) * i * t) - np.exp(tau * i * x) *
np.exp((tau2 + 4) * i * t))
l=6
N = 2 ** l
dx = 1.0 / N
dx2 = dx * dx
dt = dx2/2
t_final = 1.
x_arr = np.arange(0,1,dx)
approximate_f = np.cos(tau*x_arr)
approximate_g = -i*np.sin(tau*x_arr)
H = np.zeros([2*N,2*N], dtype=np.complex)
for k in range(N):
H[k,k] = -3*i*dt
H[k,k+N] = (2/dx2+1)*i*dt
if k==0:
H[k,N+1] = -i/dx2*dt
H[k,-1] = -i/dx2*dt
elif k==N-1:
H[N-1,N] = -i/dx2*dt
H[N-1,-2] = -i/dx2*dt
else:
H[k,k+N-1] = -i/dx2*dt
H[k,k+N+1] = -i/dx2*dt
### Bottom half
H[N :, :N] = H[:N, N:]
H[N:, N:] = H[:N, :N]
### Theta method. 0.5 -> Crank Nicolson
theta=0.5
A = np.eye(2*N)+H*theta
B = np.eye(2*N)-H*(1-theta)
### Precompute for faster computations
mat = np.linalg.inv(A)@B
t = 0
b = np.concatenate((approximate_f, approximate_g))
while t < t_final:
t += dt
b = mat@b
approximate_f = b[:N]
approximate_g = b[N:]
approximate_solution = np.concatenate((approximate_f, approximate_g))
#Calculate the actual solution
actual_f = solution_f(t,np.arange(0,1,dx))
actual_g = solution_g(t,np.arange(0,1,dx))
actual_solution = np.concatenate((actual_f, actual_g))
plt.figure(figsize=(7,5))
plt.plot(x_arr,actual_f.real,c="C0",label=r"$Re(f_\mathrm{true})$")
plt.plot(x_arr,actual_f.imag,c="C1",label=r"$Im(f_\mathrm{true})$")
plt.plot(x_arr,approximate_f.real,c="C0",ls="--",label=r"$Re(f_\mathrm{num})$")
plt.plot(x_arr,approximate_f.imag,c="C1",ls="--",label=r"$Im(f_\mathrm{num})$")
plt.legend(loc=3,fontsize=12)
plt.xlabel("x")
plt.savefig("num_approx.png",dpi=150)
