경직성(Stiff) 모델로 가변 스텝 솔버 탐색하기

다음 제품이 필요합니다.

이 예제에서는 푸코의 진자(Foucault pendulum) 모델에서 가변 스텝 솔버가 어떻게 동작하는지 보여줍니다. Simulink® 솔버 ode45, ode15s, ode23, ode23t를 테스트 케이스로 사용합니다. 이 문제의 해를 구하는 데는 경직성 미분 방정식이 사용됩니다. 방정식의 경직성에 관한 정확한 정의는 없습니다. 일부 수치적 방법은 경직성 방정식의 해를 구하는 데 사용하기에 불안정하며 경직성 문제에 대해 수치적으로 안정된 해를 구하기 위해서는 아주 작은 스텝 크기가 필요합니다. 경직성 문제는 빠르게 변하는 성분을 가질 수도 있고 느리게 변하는 성분을 가질 수도 있습니다.

경직성 문제의 한 예로 푸코 진자를 들 수 있습니다. 진자는 몇 초 이내에 진동을 완료하는(빠르게 변하는 성분) 반면, 지구는 하루에 한 번 축을 중심으로 회전을 완료합니다(느리게 변하는 성분). 진자의 진동면은 지구의 축 회전으로 인해 느리게 회전합니다. Modeling a Foucault Pendulum에서 푸코의 진자의 물리학에 관해 자세히 알아볼 수 있습니다.

이 시뮬레이션은 지구 표면에서 관찰자가 보는 것처럼 x-y 평면에서 진자 추의 위치를 계산합니다. 그런 다음, 진자의 총 에너지를 계산하여 다양한 Simulink 솔버의 성능을 평가하는 데 사용합니다.

푸코의 진자 모델

푸코의 진자는 아래와 같은 두 개의 연립미분방정식으로 나타낼 수 있습니다. 마찰과 공기 저항은 고려하지 않습니다(이로 인해 방정식이 매우 단순화됨). 이러한 방정식의 전체 전개는 푸코의 진자 예제에서 확인할 수 있습니다.

$\ddot{x} = 2\Omega \sin{\lambda} \dot{y} - \frac{g}{L} x$

$\ddot{y} = - 2\Omega\sin{\lambda} \dot{x} - \frac{g}{L} y$

$\Omega = \mbox{ Earth's angular velocity of rotation around its axis}$

$\lambda = \mbox{ geographic latitude}$

$g = \mbox{ acceleration of gravity}$

$x \mbox{ and } y = \mbox{ coordinates of the pendulum bob}$

모델 sldemo_solvers는 푸코의 진자를 나타내는 미분 방정식의 해를 구하는 데 사용됩니다. 이 예제는 86,400초 동안 푸코의 진자를 시뮬레이션합니다. 상수와 초기 조건은 모델 작업 공간에 저장됩니다.

가변 스텝 솔버

이 예제에서는 ode45, ode15s, ode23, ode23t 가변 스텝 솔버의 성능을 조사합니다.

솔버 성능 평가하기

솔버의 성능을 평가하는 데는 다양한 방법이 있습니다. 문제가 닫힌 형식의 해를 갖는다면 솔버 결과와 예상되는 이론적인 결과를 비교할 수 있습니다. 특정 솔버를 사용하여 모델을 시뮬레이션하는 데 소요되는 시간을 모니터링할 수도 있습니다.

그러나, 푸코의 진자 문제에 대해 이론적으로 정확한 해는 없습니다. 닫힌 형식 해의 근삿값이 있을 뿐입니다. 솔버 성능을 평가하는 데 닫힌 형식 해의 근삿값을 사용할 수는 없습니다(푸코의 진자 예제 참조).

총 에너지 보존

이 예제에서는 에너지 보존 법칙을 검증하여 솔버 성능을 평가합니다. 마찰이 없는 환경에서 진자의 총 에너지는 균일하게 유지되어야 합니다. 그러나 계산된 진자의 에너지는 제한된 수치 정확도로 인해 균일하지 않습니다.

이 예제에서는 모든 시간 스텝에서 정규화된 진자의 총 에너지를 계산합니다. 에너지의 상대 오차는 시뮬레이션 과정에서 총 에너지의 변화량과 같습니다. 에너지가 보존되므로 이상적으로는 상대 오차가 0이어야 합니다. 총 에너지는 잠재적인 에너지와 운동으로 발생하는 에너지의 합입니다. NormalizeEnergy 블록은 다음 방정식을 사용하여 정규화된 진자 에너지를 계산합니다.

$E = \frac{v_x^2 + v_y^2}{2} + g( L - \sqrt{L^2 - x^2 - y^2} )$

$E_{normalized}(t) = \frac{E(t)}{E(0)}$

$E(0) = \mbox{ initial total energy}$

이 플롯은 ode23과 ode23t를 사용하여 계산했을 때의 정규화된 에너지와 시간을 비교하여 보여줍니다. 이 특정 문제에서는 ode23t가 ode23보다 훨씬 더 정확합니다. ode23을 사용한 시뮬레이션에서는 정규화된 진자 에너지가 60% 넘게 감소합니다. 낮은 에너지를 갖는 진자는 진동 진폭이 더 낮습니다. 이 효과는 다음 플롯에서 확인할 수 있습니다. 여기서 ode23으로 계산한 진자의 진폭은 진동면이 회전할 때 감소합니다.

이 플롯은 경직성 솔버와 경직성이 없는 솔버 사이의 차이를 보여줍니다. 각 플롯은 하루 동안의 진자 추의 위치를 보여줍니다. 열다섯 번째 데이터 점마다 파란색 점으로 플로팅되어 있습니다. 검은색 점은 진자 추의 초기 위치를 나타내며 검은색 선은 초기 진자 진동면을 나타냅니다. 빨간색 선은 하루가 지난 뒤의 진동면을 나타냅니다. 노란색 선은 시뮬레이션 시간 동안 어떤 중간 점에서의 진동면을 보여줍니다. 진자의 진동면은 하루 동안 전체 회전이 완료되지 않습니다. 진동면이 얼마나 빠르게 회전하는지는 지리학적 위도에 따라 다릅니다(자세한 내용은 푸코의 진자 예제 참조). 왼쪽 플롯의 진자의 진폭은 감소하는 반면, 오른쪽 플롯의 진폭은 균일하게 유지됩니다. 1e-5의 동일한 상대 허용오차에 대해 경직성 솔버는 더 정확한 결과를 제공하지만 더 많은 실행 시간이 필요합니다.

다음 플롯은 Simulink 솔버의 더 세부적인 성능 연구를 보여줍니다. 상대 허용오차의 함수로서 상대 오차와 실행 시간이 어떻게 달라지는지 보기 위해 네 개의 솔버를 선택했습니다. 더 광범위한 솔버 테스트를 위해 sldemo_solvers_mcode.m 스크립트를 사용할 수 있습니다. 이 스크립트는 모델에서 C 코드를 생성하여 시뮬레이션의 속도를 높이며, 실행하는 데 몇 분 정도 소요될 수 있습니다.

Building RSim executable for model..
Time taken: 13.530949 seconds.

Running generated RSim executable..
Time taken: 197.034259 seconds.

이 예제에서 상대 허용오차 값이 1e-6보다 작은 경우 상대 오차가 많이 감소하지 않습니다. 이는 특정 모델에 따라 달라지는 수치적 솔버의 제한 사항입니다. 솔버의 상대 허용오차를 낮춘다고 해서 정확도가 반드시 개선되는 것은 아닙니다. 문제에 필요한 최소 정확도를 추정해야 하며 이에 맞는 솔버를 선택하여 시뮬레이션 비용의 균형을 맞춰야 합니다. 예를 들어, 푸코의 진자 추의 위치가 몇 센티미터 내에 있는지 알고 싶을 수 있습니다. 몇 미크론 이내에 있는 진자의 위치는 정확하게 측정할 수 없으므로 해당 위치를 계산하는 것은 필요하지 않습니다.

시뮬레이션의 수치적 결과와 동작은 솔버 설정에 따라 달라질 수 있습니다. 이는 경직성 문제에서는 중요한 사안입니다. 경직성 모델을 사용할 때는 최소한의 계산 비용으로 정확한 결과를 낼 수 있는 솔버를 선택해야 합니다. 가변 스텝 솔버의 상대 허용오차 또는 고정 스텝 솔버의 스텝 크기는 결과가 수치적으로 안정되도록 충분히 작아야 합니다. 어떤 솔버가 다른 솔버보다 더 효율적이더라도 경직성 솔버는 경직성 문제의 해를 구하는 데 더 적합합니다. 가변 스텝 솔버는 고정 스텝 솔버보다 더 견고합니다.