algorithm - 파이썬에서 암시 적 오일러로 pde 해결 - 잘못된 출력

L2-norm은 수렴을 테스트하는 데 유용한 지표가 될 수 있지만 문제가 무엇인지 설명하지 않으므로 디버깅 할 때 이상적이지 않습니다. 솔루션이 무조건 안정적이어야하지만 후진 오일러가 반드시 정답으로 수렴하지는 않습니다. 전방 오일러가 악명 높은 (방산 방지) 악명 높은 것처럼, 후방 오일러는 악명 높습니다. 솔루션을 플로팅하면이를 확인할 수 있습니다. 수치 해는 0으로 수렴합니다. 다음의 근사치에 대해서는 크랭크-니콜슨이 합리적인 후보입니다. 아래 코드에는 솔루션의 암시성을 조정할 수 있도록보다 일반적인 세타 방법이 포함되어 있습니다. theta = 0.5는 CN을 제공하고, theta = 1은 BE를 제공하며, theta = 0은 FE를 제공합니다. 내가 조정 한 다른 것들 :

나는 dt = dx 대신 dt = (dx ** 2)/2의보다 적절한 시간 간격을 선택했습니다. 후자는 CN을 사용하여 올바른 솔루션으로 수렴하지 않습니다.

사소한 참고 사항이지만 t_final이 dt의 배수임을 보장하지 않기 때문에 동시에 솔루션을 비교하지 않았습니다.

느린 속도에 대한 의견과 관련하여 공간 해상도를 높이면 시간 해상도도 증가해야합니다. dt = dx 인 경우에도 (1024 x 1024) * 1024 매트릭스 곱셈을 1024 번 수행해야합니다. 나는 이것이 내 컴퓨터에서 특히 오래 걸리는 것을 찾지 못했습니다. 나는 약간의 속도를 높이기 위해 불필요한 연결을 제거했지만 시간 단계를 dt = (dx ** 2)/2로 변경하면 불행히도 상황이 줄어들 것입니다. 속도가 걱정된다면 Numba로 컴파일을 시도 할 수 있습니다.

말한 것은 CN의 일관성으로 엄청난 성공을 거두지 못했습니다. t_final = 1에서 아무것도 얻으려면 N = 2 ^ 6을 설정해야했습니다. 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)

당신의 수학, 그러나 나는 제안을 제안 할 것입니다.

fxx 에 대한 직접 계산 사용  그리고 gxx  수치 적으로 불안정한 좋은 후보처럼 보입니다. 직관적으로 1 차 방법은 용어에서 2 차 실수를 할 것으로 예상됩니다. 개별 항의 2 차 실수는 해당 공식을 통과 한 후 2 차 도함수에서 일정한 순서의 실수로 나타납니다. 또한 단계 크기가 작아지면 2 차 공식으로 작은 반올림 실수도 놀라 울 정도로 큰 오류로 바뀝니다.

대신 이것을 f 의 4 가지 기능의 1 차 시스템으로 바꾸는 것이 좋습니다. fx g 및 gx . 그런 다음 해당 시스템에서 뒤로 오일러를 진행하십시오. 직관적으로,이 방법으로 1 차 방법은 2 차 실수를 생성하며, 2 차 실수는 1 차 실수를 생성하는 공식을 통과합니다. 이제 시작부터 수렴해야하며 반올림 오류 전파에 민감하지 않습니다.