비경직성(Nonstiff) ODE 풀기
이 페이지에는 ode45
를 사용한 비경직성 상미분 방정식 풀이에 대한 두 가지 예제가 포함되어 있습니다. MATLAB®에는 비경직성 ODE에 대한 여러 개의 솔버가 있습니다.
ode45
ode23
ode78
ode89
ode113
ode45
는 대부분의 비경직성 문제에서 가장 잘 동작합니다. 그러나, 조금 덜 엄격한 허용오차가 용인되는 문제나 반경직성이 있는 문제에는 ode23
이 권장됩니다. 마찬가지로, 엄격한 허용오차를 가지는 문제에서나 ODE 함수를 계산하는 데 시간이 많이 걸리는 경우 ode113
이 ode45
보다 더 효율적일 수 있습니다. ode78
및 ode89
는 안정성에 있어 정확도가 중요한 역할을 하는 긴 구간에서의 적분에 뛰어난 고차수 솔버입니다.
비경직성 솔버로 문제를 풀 때 시간이 오래 걸리거나 적분에 지속적으로 실패할 경우에는 문제가 경직성 문제일 수 있습니다. 자세한 내용은 경직성(Stiff) ODE 풀기 항목을 참조하십시오.
예제: 비경직성 반데르폴 방정식(Nonstiff van der Pol equation)
반데르폴 방정식은 2계 ODE입니다.
여기서 는 스칼라 파라미터입니다. 대입
를 수행하여 이 방정식을 연립 1계 ODE로 재작성합니다. 결과로 생성되는 1계 연립 ODE는 다음과 같습니다.
연립 ODE는 ODE 솔버가 사용할 수 있는 함수 파일로 코딩되어야 합니다. ODE 함수의 일반 함수 시그니처는 다음과 같습니다.
dydt = odefun(t,y)
즉, 이 함수는 계산에 t
를 사용하지 않는 경우에도 t
와 y
를 모두 입력값으로 받아야 합니다.
함수 파일 vdp1.m
은 를 사용하여 반데르폴 방정식을 코딩합니다. 변수
와
는
y(1)
과 y(2)
로 표현되며, 요소를 2개 가진 열 벡터 dydt
에는 및
에 대한 표현식이 포함되어 있습니다.
function dydt = vdp1(t,y) %VDP1 Evaluate the van der Pol ODEs for mu = 1 % % See also ODE113, ODE23, ODE45. % Jacek Kierzenka and Lawrence F. Shampine % Copyright 1984-2014 The MathWorks, Inc. dydt = [y(2); (1-y(1)^2)*y(2)-y(1)];
초기값 [2 0]
을 사용하고 시간 구간 [0 20]
에 대해 ode45
함수를 사용하여 ODE를 풉니다. 출력값은 시간 지점 t
의 열 벡터와 해 배열 y
입니다. y
의 각 행은 t
의 해당 행에 반환된 시간에 대응됩니다. y
의 첫 번째 열은 에 대응되고, 두 번째 열은
에 대응됩니다.
[t,y] = ode45(@vdp1,[0 20],[2; 0]);
t
에 대해 과
의 해를 플로팅합니다.
plot(t,y(:,1),'-o',t,y(:,2),'-o') title('Solution of van der Pol Equation (\mu = 1) using ODE45'); xlabel('Time t'); ylabel('Solution y'); legend('y_1','y_2')
vdpode
함수는 동일한 문제를 풀지만, 에 대한 사용자 지정 값을 받습니다. 반데르폴 방정식은
가 증가함에 따라 경직성(Stiff)이 됩니다. 예를 들어, 값
을 사용하는 경우
ode15s
와 같은 경직성 솔버를 사용하여 시스템을 풀어야 합니다.
예제: 비경직성 오일러(Euler) 방정식
외부 힘이 작용하지 않는 강체에 대한 오일러 방정식은 비경직성 문제용으로 제작된 ODE 솔버에 대한 표준 테스트 문제입니다.
방정식은 다음과 같습니다.
함수 파일 rigidode
는 이 1계 연립방정식을 정의하고 초기값 ,
,
에 대응하는 초기 조건 벡터
[0; 1; 1]
을 사용하여 시간 구간 [0 12]
에 대해 이 연립방정식을 풉니다. 로컬 함수 f(t,y)
는 이 연립방정식을 인코딩합니다.
rigidode
는 출력 인수 없이 ode45
를 호출하므로 솔버는 디폴트 출력 함수 odeplot
을 사용하여 각 스텝이 완료된 후 해에 해당하는 점을 자동으로 플로팅합니다.
function rigidode %RIGIDODE Euler equations of a rigid body without external forces. % A standard test problem for non-stiff solvers proposed by Krogh. The % analytical solutions are Jacobian elliptic functions, accessible in % MATLAB. The interval here is about 1.5 periods; it is that for which % solutions are plotted on p. 243 of Shampine and Gordon. % % L. F. Shampine and M. K. Gordon, Computer Solution of Ordinary % Differential Equations, W.H. Freeman & Co., 1975. % % See also ODE45, ODE23, ODE113, FUNCTION_HANDLE. % Mark W. Reichelt and Lawrence F. Shampine, 3-23-94, 4-19-94 % Copyright 1984-2014 The MathWorks, Inc. tspan = [0 12]; y0 = [0; 1; 1]; % solve the problem using ODE45 figure; ode45(@f,tspan,y0); % -------------------------------------------------------------------------- function dydt = f(t,y) dydt = [ y(2)*y(3) -y(1)*y(3) -0.51*y(1)*y(2) ];
rigidode
함수를 호출하여 비경직성 오일러 방정식을 풉니다.
rigidode title('Solution of Rigid Body w/o External Forces using ODE45') legend('y_1','y_2','y_3','Location','Best')
참고 항목
ode45
| ode23
| ode78
| ode89
| ode113