스펙트럼법으로 1차원 선형 대류 방정식 by Fortran을 실현해 보았다

개시하다
유한원법을 사용하여 구조계산을 해 왔다@soy_bot.
최근에 시간 서열 분석을 하는 작업이 매우 많아서 도구로 때때로 부립엽 분석을 사용하고 재분석을 한다.
특히 고속 부립엽 변환(Fast Fourier Transform·FFT)은 사용 빈도가 많아 학습과 함께 실장 게임도 진행한다.
이런 상황에서 보스로부터 FFT로 미분방정식을 구해내는 스펙트럼법이 있다고 들었는데, 우연히 실험실에서 텍스트가 떨어져서 실장하고 놀기로 했다.
나는 석강 선생님의 스펙트럼법 교과서에서 공부한다.
https://www.amazon.co.jp/%E3%82%B9%E3%83%9A%E3%82%AF%E3%83%88%E3%83%AB%E6%B3%95%E3%81%AB%E3%82%88%E3%82%8B%E6%95%B0%E5%80%A4%E8%A8%88%E7%AE%97%E5%85%A5%E9%96%80-%E7%9F%B3%E5%B2%A1-%E5%9C%AD%E4%B8%80/dp/4130613057
그래서 이 글의 내용은 본문의 3장 정도의 내용에 따라 쓴 것이다.
자세한 설명은 이 책을 참조하시오.
다음은 그 비망록이다.
스펙트럼법에 대해 나는 문외한이어서 무슨 잘못이 있을까 봐 두렵다.
무엇이 스펙트럼법입니까?
미분방정식의 근사해에서 우리는 기함수에서 부립엽급수와 절비요리사 다항식 등 주기함수의 총칭을 사용하는 것으로 이해한다.
또한 적당한 기함수를 사용함으로써 에너지, 동량과 위상속도를 엄격하게 보존할 수 있는 우수한 해법이다.
일차 선형 대류 방정식

\frac{\partial u(x,t) }{\partial t} + c \frac{\partial u(x,t) }{\partial x} = 0

이렇게 표시된 미분 방정식.
이번에는 아래의 초기 조건과 경계 조건을 간단하게 설정하기 위해 대류 속도 c가 간단하기 때문에 1로 설정합니다.
(지적해 주셔서 감사합니다.)

u(x,0) = f(x) \hspace{1cm} (0 \le x \le 2 \pi)

u(x,t) = u(x +2 \pi, t)

스펙트럼법에 근거한 해법
우선 미지의 함수u(x,t)를 이산 부립엽과 비슷하게 한다.

u(x,t) = \sum_{-N}^{N} \hat{u}_{k}(t) {e}^{ikx}

또 다음 공식에 따라 잔차R를 준다.

R = \frac{\partial u(x,t) }{\partial t} + \frac{\partial u(x,t) }{\partial x} 

근사 공식을 이 잔차에 대입하고 가중치w를 더해 전체 영역에서 포인트를 매기면 다음과 같은 방정식을 얻는다.

\int_{0}^{2 \pi}  \sum_{-N}^{N} w\frac{\partial \hat{u}_{k}(t)}{\partial t} {e}^{ikx} + \sum_{-N}^{N} \hat{u}_{k}(t) w \frac{\partial {e}^{ikx} dx}{\partial x} = 0

또한 아래의 방정식을 사용하여 권중을 가중시키고,

w = e^{-ikx}

가중 잔차식은 다음과 같다.

\int_{0}^{2 \pi}  \sum_{-N}^{N} \frac{\partial \hat{u}_{k}(t)}{\partial t} + \sum_{-N}^{N} ik \hat{u}_{k}(t) dx= 0

이 공식을 정리하고 삼각함수의 정교성을 사용하면 다음과 같은 공식을 얻을 수 있다.

\frac{\partial \hat{u}_k (t)}{\partial t} + ik\hat{u}_k(t) = 0 

이것은 스펙트럼 시간의 발전을 통제하는 상미분 방정식이어서 나는 매우 감동했다.
따라서 초기값을 이산 부립엽 변환하여 초기 스펙트럼을 얻은 다음에 상기 미분방정식 구해 시간에 따라 발전하면 임의의 시간에 부립엽 스펙트럼을 얻을 수 있다.
구체적으로 말하면

\hat{u}_k(0) = \frac{1}{2 \pi} \int_0^{2 \pi} f(x) e^{-ikx} dx 

계산 초기 스펙트럼,
상기 상미분 방정식을 넣고 계산함으로써 다음과 같은 방정식을 얻는다.

\hat{u}_k(t) = \hat{u}_k(0) e^{-ikt}

그리고 시간에 따라 이 부립엽 스펙트럼을 역이산 부립엽으로 변환하면 근사해를 얻을 수 있다.

u(x,t) = \sum_{k=-N}^{N} \hat{u}_k(t) e^{ikx}

또한 부립엽 스펙트럼 시간 발전의 엄격한 해결 방안에 착안하면 이 공식은 다음과 같은 공식으로 바꿀 수 있다.

u(x,t) = \sum_{k=-N}^{N} \hat{u}_k(0) e^{-ikt} e^{ikx} =\sum_{k=-N}^{N} \hat{u}_k(0)e^{ik(x-t)}

이는 이산화에 사용되는 모든 주파수의 위상 속도가 엄격해와 같다는 것을 의미한다.
그나저나 중심차분으로 비슷하면 안 되고 고주파 성분의 위상 속도가 늦어진다.
Fortran 기반 코드
그럼 구체적으로 계산해 봅시다.
직사각형 파의 초기 값이 어떻게 될지 궁금해서 해봤어요.
어쨌든
(1) 초기값을 이산 부립엽 변환하여 부립엽 스펙트럼을 구한다.
(2) 시간 진화 방정식

\hat{u}_k(t) = \hat{u}_k(0) e^{-ikt}

(인용하는 의식이 다르다. 고맙다.)
미래의 부립엽 스펙트럼을 구하다.
(3) 역이산 부립엽 변환을 통해u(x,t) 얻을 수 있다.
그러니까
다음 직사각형파의 부립엽 스펙트럼은 이번에 검증을 위해 수동으로 계산한 것이다.

f(x) = 1 \hspace{1cm} \left(\frac{4}{5}\pi \le x \le \frac{6}{5} \pi \right), 

f(x) = 0 \hspace{1cm} \left(otherwise \right)

이렇게 하면 이런 함수를 얻을 수 있다.

\hat{u}_k(0) = \frac{1}{\pi k} {e}^{-ik\pi} {\rm sin}\frac{\pi}{5}k

그리고 코드에서만 (1), (2), (3) 를 실행합니다.

program main

use plantfem
implicit none

integer(int32) ,parameter:: N = 2**9
complex(real64),allocatable :: x(:),uk0(:),uk(:),ut(:)
integer(int32) :: i,k,j
real(real64) :: T

type(Math_) :: Math
type(IO_)   :: f
type(Console_) :: console

! 初期時刻を0として，何秒後の解がほしいか．
call console%log("Input Time: T")
call console%read()
T = freal(console%line)


! (1) 離散フーリエ変換により，初期値のフーリエスペクトルを求める．
allocate(uk0(-N:N) )
do k=-N,N
    if(k==0)then
        uk0(k) = 1.0d0/Math%PI*exp(-Math%i*dble(k)*Math%PI )*Math%PI/5.0d0
        cycle  
    endif
    uk0(k) =1.0d0/Math%PI/dble(k)*exp(-Math%i*dble(k)*Math%PI )*sin(Math%PI/5.0d0*dble(k) ) 
enddo


! (2) ある時刻T後のフーリエスペクトルを求める．
allocate(uk(-N:N) )
do k=-N,N
    uk(k) = uk0(k)*exp(-Math%i*dble(k)*T )
enddo

! (3) 逆フーリエ変換(by DFT)
x = linspace([0.0d0,2.0d0*Math%PI],2*N+1)
ut = zeros(2*N+1)
do i=1,size(x)
    do k=-N,N
        ut(i) = ut(i) + uk(k)*exp(Math%i*dble(k)*x(i) )
    enddo
enddo

! 出力
call f%open("ut.txt","w")
call f%write(real(x),real(ut))
call f%close()
call f%plot(option="with lines")

end program

결과는 이렇다.




지부스 현상으로 인한 오버샷/underrshoot를 볼 수 있지만 에너지, 동량, 위상 속도가 보존되어 있다.
그러나 이번에 사용된 DFT는 매우 느려서 보통 FFT를 사용한다.
여기서 이번 파수는 2*N+1이며, 통상 FFT는 2^N으로 진행되기 때문에 약간의 노력이 필요하다.
이것은 FFT 버전을 사용했으니, 나는 즉시 투고할 것이다.
최후
이번 샘플 코드에는 포란 기능 외에 독자 확장 기능도 추가됐다.
이 기능들은planetFEM을 설치한 후에 사용할 수 있으니 반드시 사용하세요!
그럼, 즐거운 생활 되세요!