[MATLAB] 미분방정식의 수치 적분 with ode45

[MATLAB] 미분방정식의 수치 적분 with ode45

Introduction

잡담

안녕하세요~ BestEng입니다~ 오늘은 4대 역학에서 잠시 벗어나서, 오늘은 다른 주제에 대해 다뤄보겠습니다. 바로, 진화하는 공대생의 계산기의 최종 보스..(;;) MATLAB 인데요, 기본적인 인터페이스나 문법등은 알고 있다고 가정하고(그..그래도 되겠죠?) 약간의 응용에 대해 다뤄보려고 합니다.

Numerical Analysis?

공과대학 학생들이라면 공학수학 시간(참새와 함께하는 공학수학에서도..)에 지겹도록 푸는 1차, 2차 미분방정식! 여러가지 해들을 구하면서 곯머리를 많이 앓으셨을 텐데요~ 사실 이렇게 손으로 풀이가 가능할 정도로 단순한 미분 방정식들은 잘 없답니다. 특히, 비선형 항이 들어가면 더욱 난해해지죠.

하지만 공학에서는 모든 미분 방정식의 해를 직접적으로 구할 필요가 없는 경우가 많습니다. 특정한 경계 조건에서 값들이 어떻게 변화하는지만 분석하면 되는 것이죠. 예를 들어, 비행기 날개 주위의 유동을 손으로 계산한다? 불가능한 일입니다. 하지만 컴퓨터를 이용하면 우리가 설정한 경계 조건에서 시간에 따른 양력과 항력은 계산해낼 수 있습니다.

이것이 바로 수치 해석(Numerical Analysis)의 기본이며, 그 중 가장 간단한 형태의 미분방정식인 ODE(Ordinary differential equation)를 MATLAB을 이용하여 수치 적분(Numerical Integration)을 통해 해를 구하는 방법을 알아보겠습니다.

Purpose

이번 포스팅의 목표는 다음과 같습니다.

• 수치 적분의 이해
• ode45의 사용법 익히기
• 실전 예제 (참새와 함께하는 공학수학 - ODE 편 - #2-2nd order ODE(3.Problems))
• 깨알 MATLAB Tip!

그럼 시작해 보도록 하죠..!

Numerical Integration

수치 적분은 그냥 미분 방정식의 해를 수치적으로 구하는 것이구나~ 하면 됩니다. 기본적으로 수치 적분은 다음과 같은 과정을 통해 풀게 됩니다.

1. 고차 미분 방정식(High-order differential equation)을 $\dot{x}=f(x)$와 같은 꼴로 변환
2. 초기 값(or 경계 값) $x_{0}$를 위의 미분 방정식에 대입
3. 짧은 시간 간격 $\Delta t$ 동안의 $x$의 변화량 $\Delta x$를 계산
4. $x(t+\Delta t) = x(t) + \Delta x$로 업데이트
5. 3~4의 반복

여기서 $x$를 업데이트 할 때 필요한 $\Delta x$를 구하는 방법에 따라서 다양한 방법이 파생됩니다. 예를 들어, Runge-Kutta, Adams 방법 등이죠. 이번 포스팅에서는 가장 일반적으로 많이 쓰이는 MATLAB 수치 적분 함수 ode45의 사용법에 대해 알아보도록 합시다.

ode45

우선 ode45 함수의 구성부터 살펴봅시다. 아래는 ode45 함수의 MATLAB 내 도움말 중 일부입니다.

[T,Y] = solver(odefun,tspan,y0);
[T,Y] = solver(odefun,tspan,y0,options);

solver에는 우리가 쓸 ode45가 들어가고, odefun에는 우리가 풀 미분 방정식, tspan은 T의 초기 값과 마지막 값, y0에는 Y의 초기 값이 들어가면 됩니다.

Simple example

조금 더 자세히 살펴보면, odefun 자리에는 t와 y를 인수로 가지는 함수가 들어가게 됩니다. 예를 들면, 다음과 같이 쓸 수 있겠죠.

function res = myfun(t,y)
res = 2*y;

위 함수는 다음과 같은 미분 방정식을 표현한 것입니다.
$$\dot{y}=2y$$

초기 조건을 $y(0) = 1$로 두고, 시간을 $t=0 \sim1$로 두면 수치 적분은 다음과 같이 수행할 수 있습니다.

[T,Y] = ode45(myfun,[0 1],1);

혹은 익명 함수(Anonymous function)를 이용해 다음과 같이 함수 파일을 따로 저장할 필요 없이 계산할 수도 있습니다.

[T,Y] = ode45(@(t,y) 2*y,[0 1],1);

@(t,y) 2*y의 의미는 t와 y를 인자로 하고 이로부터 2*y를 계산하여 반환하는 익명 함수를 말합니다. 보다시피 따로 함수 파일을 작성할 필요 없이 명령 창에서도 작동하기 때문에 짧은 코드에서 매우 유용하게 쓸 수 있죠.

위 수치 적분 결과 값을 다음과 같이 그래프로 그려보죠.

plot(T,Y);

비록 주어진 미분 방정식이 간단한 형태라 손으로도 풀 수 있지만, 위와 같이 해를 해석적으로 구하지 않고도 시간 $t$에 따라 상태 $y$가 어떻게 변하는 지, 알 수 있답니다.

자 이제 ode45로 어떻게 수치 적분을 하는지 감이 오시나요?

실전 예제

실전 예제는 참새와 함께하는 공학수학 - ODE 편 - #2-2nd order ODE(3.Problems) 중에서 얼핏 보기에 풀이가 가장 긴 문제 2.0.2를 다뤄보도록 하겠습니다.

풀이 - 참새와 함께하는 공학수학 - ODE 편 - #2-2nd order ODE(4.문제풀이)

문제에서 주어진 2차 미분 방정식은 다음과 같네요.

$$x \frac{d^{2}y}{dx^{2}}+2\frac{dy}{dx}+xy=0$$

문제에서는 해가 하나 주어지고 나머지를 구하라고 했지만, 우리에겐 강력한 ode45가 있으니 바로 해를 구해봅시다. 우선 주어진 2차 미분 방정식을 1차 미분 방정식으로 바꿔주기 위해 새로운 상태 변수 벡터 $[y, \dot{y}]^{T}$를 잡으면, 위 식은 다음과 같이 쓸 수 있습니다.

$$\begin{Bmatrix} \dot{y} \\ \ddot{y} \end{Bmatrix} = \begin{bmatrix} 0 & 1 \\ -1 & -\frac{2}{x}\end{bmatrix} \begin{Bmatrix} y \\ \dot{y} \end{Bmatrix}$$

따라서 다음과 같은 MATLAB 코드를 구성할 수 있겠네요.

% Boundary Conditions
x0 = 1; xf = 10;

figure();
hold on;
for i = 1:10
    % Random initial values
    [X,Y] = ode45(@(x,y) [y(2);-y(1)-2/x*y(2)], [x0, xf], rand(2,1));
    plot(X,Y(:,1));
end
grid on;
xlabel('x'); ylabel('y');
hold off;

여기서, 임의의 초기 값을 생성하기 위해 for문 안에 rand 함수를 이용하여 총 10 번의 임의의 초기 값을 만들고 이를 각각 $x=1\sim10$에서 수치 적분 하였습니다. tspan(혹은 $x$)의 초기 값이 0이 아니라 1인 이유는 위 식에서 $\ddot{y}$를 $x=0$에서 계산을 할 수 없기 때문입니다.

이를 실행해 보면 임의의 초기 값에 대해 아래와 같은 $y=y(x)$ 그래프가 나옵니다.

어떄요, 참 쉽-죠?

여기서 수치적으로 적분한 결과와 참새가 열심히 손으로 푼 결과를 비교해 보도록 하겠습니다. 주어진 미분 방정식의 해석적 해는 다음과 같다고 하네요.

$$y(x) = \frac{1}{x}\left(c_{1}cos(x)+c_{2}sin(x)\right)$$

수치 해와 해석 해를 비교하기 위해 초기 값을 하나 선택하고, 둘을 같이 그려주는 MATLAB 코드를 다음과 같이 구성해 봅시다.

% Boundary Conditions
y0 = rand(2,1);
x0 = 1; xf = 10;

% Numerical Integration
[X,Y] = ode45(@(x,y) [y(2);-y(1)-2/x*y(2)], [x0, xf], y0);

% Analytic Solution
c = [cos(x0)/x0 sin(x0)/x0; (-sin(x0)-cos(x0))/x0^2 (cos(x0)-sin(x0))/x0^2]\y0;
Y_a = 1./X.*(c(1).*cos(X)+c(2).*sin(X));

figure()
% Plot Numerical and Analytic Solution
subplot(2,1,1);
hold on;
plot(X,Y(:,1),'o');
plot(X,Y_a);
grid on;
ylabel('y');
legend('Numerical Solution', 'Analytic Solution');
hold off;

% Error
subplot(2,1,2);
plot(X,Y(:,1)-Y_a);
grid on;
legend('Error')
xlabel('x'); ylabel('\Delta y');

여기서 9번 째 줄은 해석 해의 초기 값을 주어진 조건으로부터 구하는 식입니다.

이를 실행하면 다음과 같은 그래프를 그릴 수 있습니다.

보시면 수치 해와 해석 해가 거의 같은 것을 알 수 있죠?

Conclusion

이번 포스팅에서는 수치 해석 중 수치 적분에 대해 실전 예제를 곁들여서 조금 설명을 해 보았습니다. 여러 문제를 풀면서 직접 코딩도 해 보며 MATLAB에 조금씩 익숙해지면 좋을 것 같네요~

앞으로는 시간날 때마다 참새가 열심히 손으로 푼 문제들을 수치적으로는 어떻게 푸는지, 간단히 풀어보는 포스팅을 하도록 할게요~

그럼 20000~~