Вы реализовали числовой решатель для линейного (без произведения p) уравнения в частных производных, которое, кажется, работает для ряда временных шагов, а затем взрывается с быстрыми колебаниями.
Математически происходит следующее: ваша дискретная система уравнений имеет решения, которые экспоненциально растут со временем. Числовой шум в наименее значащих цифрах ваших чисел с плавающей запятой будет экспоненциально расти, пока не примет решение, которое вы ищете.
Если P - вектор (массив из 2N значений), состоящий из значения давления N на текущем и предыдущем временных шагах, все сложение и вычитание значений в вашем коде можно переписать в матричной форме (представляющей один временной шаг):
P := D @ P
где D матрица формы (2N, 2N). Если P0 было начальным состоянием системы, то состояние системы после n временных шагов будет
Pn = np.linalg.matrix_power(D, n) @ P0
Здесь matrix_power - матричное возведение в степень, т.е. matrix_power(D, 3) == D @ D @ D. Если D имеет собственные значения с модулем> 1, то будут экспоненциально растущие решения. Ваши данные показывают, что наибольшее собственное значение составляет около 1000, что означает, что оно выросло в 1e + 18 раз за 6 временных шагов, тогда как физически релевантный собственный вектор составляет около 0,9. Обратите внимание, что числа с плавающей запятой имеют точность около 1e-16.
Вы можете проверить собственные значения, используя np.linalg.eigvals, желательно с небольшой матрицей, например N=100. Если вам повезет, то уменьшения размеров шага dx и dt может быть достаточно, чтобы собственные значения оказались ниже единицы.
Быстрое и грязное решение
Это быстро реализовать, но не очень эффективно. Вы можете попытаться удалить большие собственные значения из D:
import numpy as np
evals, evecs = np.linalg.eig(D)
evals_1 = np.where(np.abs(evals) > 1, 0, evals)
D_1 = evecs @ np.diag(evals_1) @ np.linalg.inv(evecs)
Например, с небольшой матрицей (собственные значения 0,9 и 1000):
D = np.array([[500.45, -499.55],
[-499.55, 500.45]])
# D_1 calculated as above.
D_1 = np.array([[0.45, 0.45],
[0.45, 0.45]])
P0 = [1, 1+8e-16]
powm = np.linalg.matrix_power
ns = np.arange(7)
Ps_D = np.array([powm(D, n) @ P0 for n in ns])
Ps_D1 = np.array([powm(D_1, n) @ P0 for n in ns])
import matplotlib.pyplot as plt
fig, ax = plt.subplots()
ax.semilogy(ns, Ps_D[:, 1], 'ro-', label='Using D')
ax.semilogy(ns, Ps_D1[:, 1], 'bs-', label='Using D_1')
ax.set_xlabel('Iteration number')
ax.set_ylabel('P[1]')
ax.grid()
ax.legend()
fig.show()
If you went through the trouble of finding the eigenvalues and eigenvectors (which can be slow for large matrices), you can speed up the subsequent calculations quite a bit by replacing matpow(D_1, n) by
Pn = evecs @ np.diag(evals_1**n) @ np.linalg.inv(evecs) @ P0
Quick sanity check for tuning step sizes
It is possible to judge whether it is likely that there are eigenvalues greater than one by setting your initial state to
P = np.zeros(N)
P[N//4] = 1
and executing one single time step. If this results in values >1 in P, you are likely to have too large eigenvalues. (The new value for P is actually the column D[:, N//4] of the matrix).
Implicit methods
Rather than bluntly zeroing large eigenvalues as described above, it may be better to resort to an неявный метод . Для этого необходимо записать матрицу D. Если P - текущее состояние, а Pnew - состояние для следующего временного шага, тогда вы начинаете с уравнения
Pnew = D @ (0.5*(P + Pnew))
Решите это матричное уравнение для Pnew:
# Pnew = inv(1 - 0.5 D) (0.5 D) P
nD = D.shape[0]
D_2 = np.linalg.inv(np.eye(nD) - 0.5*D) @ (0.5*D)
# time step:
Pnew = D_2 @ P
Это метод Крэнка-Николсона . Это требует, чтобы вы построили матрицу D. Применение к тому же примеру, что и выше, дает:
D_2 = np.array([[-0.09191109, 0.91009291],
[ 0.91009291, -0.09191109]])
evals = np.array[ 0.81818182, -1.00200401]
Это уменьшает наибольшее собственное значение с 1000 до -1,002, что лучше, но все же больше 1 (по абсолютной величине), поэтому, начиная с ошибки округления на 1e-15 он обгонит желаемое решение примерно через 34k временных шагов. Более того, это уменьшает собственное значение 0,9 до 0,8. Вам нужно будет выяснить, какой из них лучше описывает физику. К сожалению, численное решение уравнений в частных производных сложно.
РЕДАКТИРОВАТЬ: Лутц Леманн указал, что работа с инверсиями больших матриц (результат точной сетки по x) не является хорошей идеей. Обратите внимание, что D - это разреженная матрица, лучше представленная как scipy.sparse.csr_matrix, csc_matrix или dia_matrix, чем np.array. Затем вы должны решить
(1 - 0.5*D) @ P_new = 0.5*D @ P
, которое представляет собой матричное уравнение в форме Ax = b с векторами x и b:
A = scipy.sparse.identity(nD) - 0.5*D
for every time step:
b = D @ (0.5*P)
P = scipy.sparse.linalg.spsolve(A, b)
Обратите внимание, что spsolve все еще может быть медленный. Это может помочь расположить P с чередованием текущего и предыдущего значений давления, чтобы D было полосно-диагональным и могло быть сохранено как разреженное dia_matrix.
Изменить:
Следовательно, наиболее эффективный подход будет таким:
import scipy.sparse as sp
n = 200 # Example - number of x points
nD = 2*n # matrix size
D = sp.dok_matrix((nD, nD))
# Initialize pressures
P = np.zeros(nD)
# ... initialize P and nonzero matrix elements here ...
# Prepare matrices
D = sp.csc_matrix(D) # convert to efficient storage format
A = sp.identity(nD) - 0.5*D
A_lu = sp.linalg.splu(A)
for _i in range(num_time_steps):
b = D @ (0.5*P)
P = A_lu.solve(b)