홈>
현재 진행중인 상황과 문제점을 정확하게 설명하고 설명하겠습니다.
이것은 약간 교묘하고 그래서 라텍스를 지원하지 않으므로 슬프게도 이미지에 의지해야했습니다. 괜찮 길 바랍니다.
반전 된 이유를 모르겠습니다. 죄송합니다. 여하튼, 이것은 우리가 A와 b를 아는 선형 시스템 Ax = b이므로 x를 찾을 수 있습니다. 이는 다음 시간 단계에서의 근사치입니다. 시간 t_final까지 계속합니다.
이것은 코드입니다
import numpy as np
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))
for l in range(2, 12):
N = 2 ** l #number of grid points
dx = 1.0 / N #space between grid points
dx2 = dx * dx
dt = dx #time step
t_final = 1
approximate_f = np.zeros((N, 1), dtype = np.complex)
approximate_g = np.zeros((N, 1), dtype = np.complex)
#Insert initial conditions
for k in range(N):
approximate_f[k, 0] = np.cos(tau * k * dx)
approximate_g[k, 0] = -i * np.sin(tau * k * dx)
#Create coefficient matrix
A = np.zeros((2 * N, 2 * N), dtype = np.complex)
#First row is special
A[0, 0] = 1 -3*i*dt
A[0, N] = ((2 * dt / dx2) + dt) * i
A[0, N + 1] = (-dt / dx2) * i
A[0, -1] = (-dt / dx2) * i
#Last row is special
A[N - 1, N - 1] = 1 - (3 * dt) * i
A[N - 1, N] = (-dt / dx2) * i
A[N - 1, -2] = (-dt / dx2) * i
A[N - 1, -1] = ((2 * dt / dx2) + dt) * i
#middle
for k in range(1, N - 1):
A[k, k] = 1 - (3 * dt) * i
A[k, k + N - 1] = (-dt / dx2) * i
A[k, k + N] = ((2 * dt / dx2) + dt) * i
A[k, k + N + 1] = (-dt / dx2) * i
#Bottom half
A[N :, :N] = A[:N, N:]
A[N:, N:] = A[:N, :N]
Ainv = np.linalg.inv(A)
#Advance through time
time = 0
while time < t_final:
b = np.concatenate((approximate_f, approximate_g), axis = 0)
x = np.dot(Ainv, b) #Solve Ax = b
approximate_f = x[:N]
approximate_g = x[N:]
time += dt
approximate_solution = np.concatenate((approximate_f, approximate_g), axis=0)
#Calculate the actual solution
actual_f = np.zeros((N, 1), dtype = np.complex)
actual_g = np.zeros((N, 1), dtype = np.complex)
for k in range(N):
actual_f[k, 0] = solution_f(t_final, k * dx)
actual_g[k, 0] = solution_g(t_final, k * dx)
actual_solution = np.concatenate((actual_f, actual_g), axis = 0)
print(np.sqrt(dx) * np.linalg.norm(actual_solution - approximate_solution))
작동하지 않습니다. 적어도 처음에는 아니고, 천천히 시작해서는 안됩니다. 무조건 안정적이어야하고 정답에 수렴해야합니다.
여기서 무엇이 잘못 되었나요?
- 답변 # 1
- 답변 # 2
당신의 수학, 그러나 나는 제안을 제안 할 것입니다.
fxx
에 대한 직접 계산 사용 그리고gxx
수치 적으로 불안정한 좋은 후보처럼 보입니다. 직관적으로 1 차 방법은 용어에서 2 차 실수를 할 것으로 예상됩니다. 개별 항의 2 차 실수는 해당 공식을 통과 한 후 2 차 도함수에서 일정한 순서의 실수로 나타납니다. 또한 단계 크기가 작아지면 2 차 공식으로 작은 반올림 실수도 놀라 울 정도로 큰 오류로 바뀝니다.대신 이것을
f
의 4 가지 기능의 1 차 시스템으로 바꾸는 것이 좋습니다.fx
g
및gx
. 그런 다음 해당 시스템에서 뒤로 오일러를 진행하십시오. 직관적으로,이 방법으로 1 차 방법은 2 차 실수를 생성하며, 2 차 실수는 1 차 실수를 생성하는 공식을 통과합니다. 이제 시작부터 수렴해야하며 반올림 오류 전파에 민감하지 않습니다.
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 또는 기타 선형 다단계 방법을 구현하여이를 개선 할 수 있습니다.
플롯이있는 코드는 다음과 같습니다.