미분대수 방정식(DAE) 풀기

미분대수 방정식이란?

미분대수 방정식은 종속 변수에 대한 하나 이상의 도함수가 방정식에 없는 미분 방정식 유형입니다. 도함수 없이 방정식에 나타나는 변수를 대수(Algebraic)라고 하며, 대수 변수가 존재한다는 것은 양함수 형식 $y' = f (t, y)$ 로 방정식을 작성할 수 없음을 의미합니다. 그 대신, 다음 형식의 DAE를 풀 수 있습니다.

ode15s 솔버와 ode23t 솔버는 다음 형식의 반명시적(semi-explicit) DAE를 포함하여 특이 질량 행렬 $M (t, y) y' = f (t, y)$ 를 가지는 지수-1 선형 음함수 문제를 풀 수 있습니다.
$\begin{matrix} y' = f (t, y, z) \\ 0 = g (t, y, z) . \end{matrix}$
이 형식에서 대수 변수가 존재하면 주대각선에 0이 하나 이상 있으므로 특이 질량 행렬이 생성됩니다.
$M y' = (\begin{matrix} y'_{1} & 0 & \dots & 0 \\ 0 & y'_{2} & 0 & ⋮ \\ ⋮ & 0 & ⋱ & 0 \\ 0 & \dots & 0 & 0 \end{matrix}) .$
기본적으로, 솔버는 질량 행렬의 특이점을 자동으로 테스트하여 연립 DAE를 감지합니다. 사전에 특이점을 알고 있는 경우 odeset의 MassSingular 옵션을 'yes'로 설정할 수 있습니다. DAE에서는 odeset의 InitialSlope 속성을 사용하여 솔버에 $y'_{0}$ 에 대한 초기 조건의 추측값을 제공할 수도 있습니다. 솔버에 대한 호출에서 $y_{0}$ 에 대해 일반적인 초기 조건을 지정하는 작업과 더불어 이 작업을 수행합니다.
ode15i 솔버는 완전 음함수 형식의 더욱 일반적인 DAE를 풀 수 있습니다.
$f (t, y, y') = 0 .$
완전 음함수 형식에서 대수 변수가 존재하면 특이 야코비 행렬이 생성됩니다. 이는 해당 변수의 도함수가 방정식에서 나타나지 않으므로 행렬에서 하나 이상의 열에 포함된 모든 요소가 0이라는 것이 보장되기 때문입니다.
$J = \partial f / \partial y' = (\begin{matrix} \frac{\partial f_{1}}{\partial y'_{1}} & \dots & \frac{\partial f_{1}}{\partial y'_{n}} \\ ⋮ & ⋱ & ⋮ \\ \frac{\partial f_{m}}{\partial y'_{1}} & \dots & \frac{\partial f_{m}}{\partial y'_{n}} \end{matrix})$
ode15i 솔버를 사용할 경우 $y'_{0}$ 과 $y_{0}$ 모두에 대해 초기 조건을 지정해야 합니다. 또한, 다른 ODE 솔버와 달리 ode15i에는 추가 입력값을 받도록 방정식을 인코딩하는 함수 odefun(t,y,yp)가 필요합니다.

물리학의 보존 법칙이 종종 $x + y + z = 0$ 과 같은 형식을 가지기 때문에 DAE는 매우 다양한 시스템에서 나타납니다. x, x', y, y'가 방정식에서 명시적으로 정의된 경우, 이러한 보존 방정식을 z'에 대한 표현식 없이 z에 대해 충분히 풀 수 있습니다.

모순 없는 초기 조건

DAE를 푸는 경우 $y'_{0}$ 과 $y_{0}$ 모두에 대해 초기 조건을 지정할 수 있습니다. ode15i 솔버를 사용할 경우 두 초기 조건이 입력 인수로 지정되어야 합니다. ode15s 솔버와 ode23t 솔버를 사용할 경우 $y'_{0}$ 의 초기 조건은 선택 사항입니다(단, odeset의 InitialSlope 옵션을 사용하여 지정할 수 있음). 두 경우 모두, 지정하는 초기 조건이 해를 구하려는 방정식과 일치하지 않을 수 있습니다. 다른 초기 조건과 서로 충돌하는 초기 조건을 일관되지 않은 초기 조건이라고 합니다. 초기 조건 처리는 솔버에 따라 다릅니다.

ode15s와 ode23t — $y'_{0}$ 에 대한 초기 조건을 지정하지 않으면 솔버는 사용자가 $y_{0}$ 에 대해 제공하는 초기 조건을 기반으로 모순 없는 초기 조건을 자동으로 계산합니다. $y'_{0}$ 에 대해 모순 있는 초기 조건을 지정하면 솔버가 해당 값을 추측값으로 처리하고 이 추측값에 가까운 모순 없는 값을 계산하려고 시도한 후 문제 풀이를 계속합니다.
ode15i — 솔버에 제공하는 초기 조건이 일관되어야 하며 ode15i는 제공된 값에 대해 일관성을 검사하지 않습니다. 헬퍼 함수 decic가 이를 위해 모순 없는 초기 조건을 계산합니다.

미분 지수(Differential Index)

DAE는 특이점의 척도인 미분 지수(Differential Index)로 특징지을 수 있습니다. 방정식을 미분하여 대수 변수를 제거할 수 있으며, 충분히 여러 번 이 작업을 수행하면 방정식이 양함수 연립 ODE 형식을 취합니다. 연립 DAE의 미분 지수는 이 연립방정식을 이에 상응하는 양함수 연립 ODE로 표현하는 데 필요한 미분 횟수를 나타냅니다. 따라서, ODE의 미분 지수는 0입니다.

지수-1 DAE의 예는 다음과 같습니다.

$y (t) = k (t) .$

이 방정식의 경우, 한 번 미분하여 양함수 ODE 형식을 얻을 수 있습니다.

$y' = k' (t) .$

지수-2 DAE의 예는 다음과 같습니다.

$\begin{array}{l} y'_{1} = y_{2} \\ 0 = k (t) - y_{1} . \end{array}$

이러한 방정식의 경우, 양함수 ODE 형식으로 바꾸려면 2개의 도함수가 필요합니다.

$\begin{array}{l} y'_{1} = k' (t) \\ y'_{2} = k'' (t) . \end{array}$

ode15s 솔버와 ode23t 솔버는 지수 1의 DAE만 풉니다. 지수가 2 이상인 방정식의 경우에는 이에 상응하는 지수-1 연립 DAE로 바꾸어야 합니다. 미분하면 연립 DAE를 이에 상응하는 지수-1 연립 DAE로 항상 바꿀 수 있습니다. 대수 방정식을 해당 도함수로 바꾸면 일부 제약 조건이 제거될 수 있다는 것에 유의하십시오. 방정식에 원래 제약 조건이 더 이상 포함되지 않을 경우 수치 해가 변동될 수 있습니다.

Symbolic Math Toolbox™를 사용하는 경우 자세한 내용은 Solve Differential Algebraic Equations (DAEs) (Symbolic Math Toolbox) 항목을 참조하십시오.

비음 제약 조건 적용하기

대부분의 odeset의 옵션은 DAE 솔버 ode15s, ode23t 및 ode15i에서 예상대로 작동합니다. 그러나 한 가지 중요한 예외는 NonNegative 옵션을 사용할 때입니다. 음함수 솔버(ode15s, ode23t, ode23tb)로 질량 행렬을 갖는 문제를 풀 때에는 이러한 솔버에서 NonNegative 옵션이 지원되지 않습니다. 따라서, 특이 질량 행렬을 반드시 갖는 DAE 문제에 비음 제약 조건을 적용하기 위해 이 옵션을 사용할 수 없습니다. 자세한 내용은 [1] 항목을 참조하십시오.

로버트슨 문제(Robertson Problem)를 반명시적(Semi-Explicit) 미분대수방정식(DAE)으로 풀기

스크립트 열기

이 예제에서는 연립 ODE를 연립 미분대수 방정식(DAE)으로 다시 만듭니다. hb1ode.m에서 확인할 수 있는 로버트슨 문제는 경직성(Stiff) ODE를 푸는 프로그램을 위한 기본 테스트 문제입니다. 연립방정식은 다음과 같습니다.

$\begin{array}{cl} y'_1 &= -0.04y_1 + 10^4 y_2y_3\\ y'_2 &= 0.04y_1 -
10^4 y_2y_3- (3 \times 10^7)y_2^2\\ y'_3 &= (3 \times
10^7)y_2^2.\end{array}$

hb1ode는 정상 상태에 대해 초기 조건 , , 을 사용하여 이 연립 ODE를 풉니다. 하지만, 이 방정식은 다음과 같은 선형 보존 법칙도 충족합니다.

해 및 초기 조건에 대해 다음과 같은 보존 법칙이 성립합니다.

이 연립방정식은 보존 법칙을 사용하여 의 상태를 결정하는 방식으로 연립 DAE로 재작성할 수 있습니다. 그러면 다음과 같은 연립 DAE가 됩니다.

$\begin{array}{cl} y'_1 &= -0.04y_1 + 10^4 y_2y_3\\ y'_2 &= 0.04y_1 -
10^4 y_2y_3-(3 \times 10^7)y_2^2\\ 0 &= y_1 + y_2 + y_3 - 1.\end{array}$

이 연립방정식을 연립 ODE로 만들려면 의 1계 도함수만 필요하므로 이 연립방정식의 미분 지수는 1입니다. 따라서, 시스템을 풀기 전에 더 변환할 필요가 없습니다.

함수 robertsdae는 이 연립 DAE를 담고 있습니다. 예제를 실행하려면 robertsdae.m을 현재 폴더에 저장하십시오.

function out = robertsdae(t,y)
out = [-0.04*y(1) + 1e4*y(2).*y(3)
   0.04*y(1) - 1e4*y(2).*y(3) - 3e7*y(2).^2
   y(1) + y(2) + y(3) - 1 ];

이러한 로버트슨 문제의 수식화에 대한 전체 예제 코드는 hb1dae.m에서 확인할 수 있습니다.

ode15s를 사용하여 연립 DAE를 풉니다. y0에 대한 모순 없는 초기 조건은 명백히 보존 법칙을 기반으로 합니다. 다음과 같이 odeset을 사용하여 옵션을 설정합니다.

상수 질량 행렬을 사용하여 연립방정식의 좌변을 나타냅니다.

$\left( \begin{array}{c} y'_1\\ y'_2\\ 0 \end{array} \right) = M y'
\rightarrow M = \left( \begin{array}{ccc} 1 & 0 & 0\\ 0 & 1 & 0\\ 0 & 0 &
0 \end{array} \right)$

상대 허용오차를 1e-4로 설정합니다.
두 번째 해 성분에 대해서는 스케일이 다른 성분에 비해 급격히 달라지므로 절대 허용오차 1e-10을 사용합니다.
DAE의 자동 감지를 테스트하기 위해 'MassSingular' 옵션을 디폴트 값 'maybe'로 그대로 둡니다.

y0 = [1; 0; 0];
tspan = [0 4*logspace(-6,6)];
M = [1 0 0; 0 1 0; 0 0 0];
options = odeset('Mass',M,'RelTol',1e-4,'AbsTol',[1e-6 1e-10 1e-6]);
[t,y] = ode15s(@robertsdae,tspan,y0,options);

해를 플로팅합니다.

y(:,2) = 1e4*y(:,2);
semilogx(t,y);
ylabel('1e4 * y(:,2)');
title('Robertson DAE problem with a Conservation Law, solved by ODE15S');

참고 문헌

[1] Shampine, L.F., S. Thompson, J.A. Kierzenka, and G.D. Byrne. “Non-Negative Solutions of ODEs.” Applied Mathematics and Computation 170, no. 1 (November 2005): 556–569. https://doi.org/10.1016/j.amc.2004.12.011.

참고 항목

ode15s | ode23t | ode15i | odeset

외부 웹사이트

Solving Index-1 DAEs in MATLAB and Simulink