>

현재 진행중인 상황과 문제점을 정확하게 설명하고 설명하겠습니다.

이것은 약간 교묘하고 그래서 라텍스를 지원하지 않으므로 슬프게도 이미지에 의지해야했습니다. 괜찮 길 바랍니다.

반전 된 이유를 모르겠습니다. 죄송합니다. 여하튼, 이것은 우리가 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

    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)
    
    

  • 답변 # 2

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

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

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

  • 이전 windows subsystem for linux - Nginx 114 및 php72-fpm은 PHP를 응용 프로그램/옥텟 스트림으로 다운로드하여 반환합니다
  • 다음 graph - 너비 우선 검색 대신 JAVA 프로그램에서 깊이 우선 검색 사용