Laplace Transform Matlab으로 해보기
귀찮아서 계속 미루다가 업로드가 늦어졌다.
이번에는 라플라스 변환을 매트랩을 통해서 실습을 해보는 내용이다.
그 전에, 라플라스 변환이 공학적으로 어떤 의미를 갖는지 간단하게 생각해보고 시작하자.
라플라스 변환, 그리고 그 변환의 특수한 경우인 푸리에 변환은 시간축을 주파수축으로 바꿔주는 축변환이다.
우리는 시간이라는 domain 안에서 살고 있어서 시간축의 정보가 더 직관적이고 익숙하게 느껴지지만, 위의 그림과 같이 주파수축으로 관찰을 하였을 때 시간축에서는 얻기 힘든 정보를 얻을 수 있는 경우가 생길수도 있다.
대개 크게 말하면 푸리에 변환은 신호 및 시스템, 통신 쪽에서 많이 쓰이고 라플라스 변환은 기계공학쪽에서 많이 쓰이는 것 같다. 특히, 라플라스 변환은 직접 풀기 힘든 꼴의 미분방정식을 풀 때도 많이 쓰이고, 이때문에 미분적분학 수업에서 라플라스 변환을 다루는 경우가 많다.
라플라스 변환을 이용해서 풀어줄 때는, 주어진 미분방정식을 라플라스 변환을 통해 s-domain으로 바꿔주고, s-domain에서 정리를 해준 후에 inverse laplace transform, 라플라스 역변환을 통해 원래의 시간축으로 바꿔서 미방을 해를 구해주게 된다. 이 이후의 부분에서 매트랩에서의 laplace transform 명령어를 알아보고, 직접 이를 이용해서 미분방정식까지 풀어볼 예정이다.
syms x y
f = 1/sqrt(x);
Lf = laplace(f)
ilaplace(Lf)
가장 먼저 1/sqrt(x)라는 함수에 대해서 기본적인 명령어를 확인해보면, 위와 같은 코드로 laplace()와 ilaplace()라는 명령어로 라플라스 변환과 역변환을 해줄 수 있다는 것을 알 수 있다. 여기서 보면, 역변환 결과가 x가 아닌 t에 대한 함수로 출력이 되는 것을 확인할 수 있는데, 이는 기본적으로 원래 domain을 시간축으로 디폴트로 봐주기 때문이다. 그리고 이 부분 역시 원하는 변수로 지정을 해줄 수 있다.
다음으로 독립변수와 변환변수를 직접 지정해주는 방법이다.
syms a t
f = exp(-a*t);
Lf = laplace(f)
Lf_1 = laplace(f, y)
Lf_2 = laplace(f, a, y)
ilaplace(Lf)
ilaplace(Lf, x)
ilaplace(Lf, a, x)
각각 Lf_1은 변환변수를 y로, Lf_2는 변환변수는 y, 독립변수는 a로 지정해주었고, 역변환 역시 변환변수를 지정해준 경우와, 변환변수와 독립변수를 모두 지정해준 경우에 대한 예시이다.
다음은, 특수한 함수인 delta function과 unit step function의 laplace transform을 실행시킨 결과인데, 일단 이 두 함수를 먼저 정의하고 가도록 하겠다.
delta function의 경우는 x축에 대해서 x=0일때만 infinite의 값을 가지는 함수로 정의가 된다. 이 함수의 특징은 전체 구간에서의 적분 결과가 1로 나온다는 것이다.
unit step function의 경우는 0점을 기준으로 하여 그 이후로는 1로 constant한 값을 갖는 계단함수이다. 이 함수의 경우는 앞서서 정의한 delta function을 통해 위와 같이 정의를 해줄 수 있다.
syms t s
Ld = laplace(dirac(t), t, s)
Lu = laplace(heaviside(t), t, s)
delta function은 dirac-delta function으로도 불리며 matlab 안에서는 dirac()이라는 명령어로 저장이 되어 있다. unit step function 역시 heaviside function이라고도 하며, 매트랩에서는 heaviside()라는 명령어로 실행시킬 수 있다. 재밌는 점은, 다를 때는 모두 0이고 시간이 0일때만 무한대의 값을 갖는 delta function을 주파수축으로 변환을 시켜주면 모든 값에서 1이라는 constant의 값을 갖는 함수가 된다는 점 정도가 될 것 같다.
마지막으로, 이제 laplace transform을 통해 직접 미분방정식의 해를 구해보도록 하자.
이를 위한 예제로는, 매우 간단하게 모델링이 된 mass spring damper system에 대해서 풀어주었고, input은 x축으로 1만큼 shift된 delta function으로 주었다. 우리는 위의 그림에서의 이계미분방정식의 해를 정확히 구해주지는 않았고, 다만 이를 라플라스 변환을 시켜주어서 s-domain에서 정리를 해주었다.
syms s t
Y_s = exp(-s)/(s^2+3*s+2);
ff = ilaplace(Y_s)
f = dsolve('D2y + 3*Dy + 2*y == dirac(t-1)','y(0)==0','Dy(0)==0');
ezplot(ff, [0,10])
hold on
s = ezplot(f, [0,10])
set(s, 'linestyle', '--', 'linewidth', 2);
코드는 우리가 쉽게 laplace 변환을 통해 s-domain에서 식을 정리해줄 수 있기 떄문에 이를 역변환하여 시간축에서의 미분 방정식의 해를 구해주고, 이를 matlab 내장 함수인 dsolve()명령어를 통해 직접 미분 방정식을 풀어서 구해준 결과와 비교해주었다.
그 결과, 위의 그래프처럼 두 결과가 같게 나온다는 것을 확인할 수 있다. (다시 코드를 돌리긴 귀찮아서 원격강의 녹화한 것에서 캡쳐해서 윗부분이 살짝 잘렸다.)
이렇게 매우 기초적인 내용에 대해서 laplace transform을 matlab을 통해 직접 실습을 해 봤다.
<추가된 내용>
이 수업들을 원격으로 해서 내 유튜브 계정에 수업 영상들이 남아있었는데, 혹시나 해서 일부공개는 그대로 넣어놓고 재생목록을 만들어두었으니 도움이 될지는 모르겠지만 관심 있으면 아래에서 영상도 확인할 수 있습니다.
https://youtube.com/playlist?list=PLT8Ck3FhPFERRvX1Y1_Ukq97u6z6pkyCG