상미분 방정식(ODE) 피팅하기

라이브 스크립트 열기

이 예제에서는 두 가지 방식으로 ODE의 파라미터를 데이터에 피팅하는 방법을 보여줍니다. 첫 번째 방법에서는 초기 파라미터에 민감하게 의존하는 유명한 ODE인 로렌츠 연립방정식의 해 중 일부에 대한 등속 원형 경로의 직관적인 피팅을 보여줍니다. 두 번째 방법에서는 등속 원형 경로를 피팅하기 위해 로렌츠 방정식의 파라미터를 수정하는 방법을 보여줍니다. 응용 사례에 적합한 접근법을 모델로 삼아 미분 방정식을 데이터에 피팅할 수 있습니다.

로렌츠 방정식: 정의 및 수치 해

로렌츠 연립방정식은 연립상미분방정식입니다(로렌츠 방정식 참조). 실수 상수 $σ, ρ, β$ 에 대해 연립방정식은 다음과 같습니다.

$\begin{array}{l} \frac{d x}{d t} = σ (y - x) \\ \frac{d y}{d t} = x (ρ - z) - y \\ \frac{d z}{d t} = x y - β z . \end{array}$

어떤 민감한 연립방정식에 대한 파라미터의 로렌츠 값이 $σ = 10, β = 8 / 3, ρ = 28$ 입니다. [x(0),y(0),z(0)] = [10,20,10]에서 시작해 시간 0부터 100까지 연립방정식의 전개를 살펴봅니다.

sigma = 10;
beta = 8/3;
rho = 28;
f = @(t,a) [-sigma*a(1) + sigma*a(2); rho*a(1) - a(2) - a(1)*a(3); -beta*a(3) + a(1)*a(2)];
xt0 = [10,20,10];
[tspan,a] = ode45(f,[0 100],xt0);     % Runge-Kutta 4th/5th order ODE solver
figure
plot3(a(:,1),a(:,2),a(:,3))
view([-10.0 -2.0])

Figure contains an axes object. The axes object contains an object of type line.

전개 양상이 매우 복잡합니다. 하지만 짧은 시간 간격에서는 균일한 원을 그리는 것처럼 보입니다. 시간 간격 [0,1/10]에 대해 해를 플로팅합니다.

[tspan,a] = ode45(f,[0 1/10],xt0);     % Runge-Kutta 4th/5th order ODE solver
figure
plot3(a(:,1),a(:,2),a(:,3))
view([-30 -70])

Figure contains an axes object. The axes object contains an object of type line.

ODE 해에 원형 경로 피팅하기

원형 경로의 방정식에는 다음과 같이 여러 개의 파라미터가 있습니다.

x-y 평면에서 경로의 각도 $θ (1)$
x축을 따라 기울어진 평면의 각도 $θ (2)$
반지름 R
속도 V
시간 0으로부터의 이동 t0
공간 delta에서의 3차원 이동

이러한 파라미터의 관점에서 시간 xdata에 대한 원형 경로의 위치를 결정하십시오.

type fitlorenzfn

function f = fitlorenzfn(x,xdata)

theta = x(1:2);
R = x(3);
V = x(4);
t0 = x(5);
delta = x(6:8);
f = zeros(length(xdata),3);
f(:,3) = R*sin(theta(1))*sin(V*(xdata - t0)) + delta(3);
f(:,1) = R*cos(V*(xdata - t0))*cos(theta(2)) ...
    - R*sin(V*(xdata - t0))*cos(theta(1))*sin(theta(2)) + delta(1);
f(:,2) = R*sin(V*(xdata - t0))*cos(theta(1))*cos(theta(2)) ...
    - R*cos(V*(xdata - t0))*sin(theta(2)) + delta(2);

ODE 해에 지정된 시간에 로렌츠 연립방정식에 대한 최적 피팅의 원형 경로를 찾으려면 lsqcurvefit을 사용하십시오. 다양한 파라미터들이 적절한 한계 내에 있도록 파라미터의 범위를 정하십시오.

lb = [-pi/2,-pi,5,-15,-pi,-40,-40,-40];
ub = [pi/2,pi,60,15,pi,40,40,40];
theta0 = [0;0];
R0 = 20;
V0 = 1;
t0 = 0;
delta0 = zeros(3,1);
x0 = [theta0;R0;V0;t0;delta0];
[xbest,resnorm,residual] = lsqcurvefit(@fitlorenzfn,x0,tspan,a,lb,ub);

Local minimum possible.

lsqcurvefit stopped because the final change in the sum of squares relative to 
its initial value is less than the value of the function tolerance.

<stopping criteria details>

ODE 해의 시간에 최적 피팅의 원형 점을 로렌츠 연립방정식의 해와 함께 플로팅합니다.

soln = a + residual;
hold on
plot3(soln(:,1),soln(:,2),soln(:,3),'r')
legend('ODE Solution','Circular Arc')
hold off

Figure contains an axes object. The axes object contains 2 objects of type line. These objects represent ODE Solution, Circular Arc.

figure
plot3(a(:,1),a(:,2),a(:,3),'b.','MarkerSize',10)
hold on
plot3(soln(:,1),soln(:,2),soln(:,3),'rx','MarkerSize',10)
legend('ODE Solution','Circular Arc')
hold off

Figure contains an axes object. The axes object contains 2 objects of type line. One or more of the lines displays its values using only markers These objects represent ODE Solution, Circular Arc.

원호에 ODE 피팅하기

이제 원호에 최적으로 피팅되도록 파라미터 $σ, β, a n d ρ$ 를 수정합니다. 더 나은 피팅을 위해 초기점 [10,20,10]도 변경할 수 있도록 합니다.

그렇게 하려면 ODE 피팅의 파라미터를 받아 시간 t에서 궤적을 계산하는 함수 파일 paramfun을 작성하십시오.

type paramfun

function pos = paramfun(x,tspan)

sigma = x(1);
beta = x(2);
rho = x(3);
xt0 = x(4:6);
f = @(t,a) [-sigma*a(1) + sigma*a(2); rho*a(1) - a(2) - a(1)*a(3); -beta*a(3) + a(1)*a(2)];
[~,pos] = ode45(f,tspan,xt0);

최적의 파라미터를 구하려면 lsqcurvefit을 사용하여 새로 계산된 ODE 궤적과 원호 soln 간의 차이를 최소화하십시오.

xt0 = zeros(1,6);
xt0(1) = sigma;
xt0(2) = beta;
xt0(3) = rho;
xt0(4:6) = soln(1,:);
[pbest,presnorm,presidual,exitflag,output] = lsqcurvefit(@paramfun,xt0,tspan,soln);

Local minimum possible.

lsqcurvefit stopped because the final change in the sum of squares relative to 
its initial value is less than the value of the function tolerance.

<stopping criteria details>

이 최적화로 인해 파라미터가 얼마나 변경되었는지 확인합니다.

fprintf('New parameters: %f, %f, %f',pbest(1:3))

New parameters: 9.132446, 2.854998, 27.937986

fprintf('Original parameters: %f, %f, %f',[sigma,beta,rho])

Original parameters: 10.000000, 2.666667, 28.000000

파라미터 sigma와 beta가 약 10% 변경되었습니다.

수정된 해를 플로팅합니다.

figure
hold on
odesl = presidual + soln;
plot3(odesl(:,1),odesl(:,2),odesl(:,3),'b')
plot3(soln(:,1),soln(:,2),soln(:,3),'r')
legend('ODE Solution','Circular Arc')
view([-30 -70])
hold off

Figure contains an axes object. The axes object contains 2 objects of type line. These objects represent ODE Solution, Circular Arc.

ODE 피팅의 문제

시뮬레이션 또는 상미분 방정식 최적화하기에서 설명한 바와 같이, 최적화 함수는 수치 ODE 해의 고유한 잡음으로 인해 문제가 있을 수 있습니다. 종료 메시지나 종료 플래그가 잠재적 부정확성을 나타낼 수 있으므로 해가 이상적이지 않다는 의심이 들면 유한 차분을 변경해 보십시오. 이 예제에서 더 큰 유한 차분 스텝 크기와 중심 유한 차분을 사용해 봅니다.

options = optimoptions('lsqcurvefit','FiniteDifferenceStepSize',1e-4,...
    'FiniteDifferenceType','central');
[pbest2,presnorm2,presidual2,exitflag2,output2] = ...
    lsqcurvefit(@paramfun,xt0,tspan,soln,[],[],options);

Local minimum possible.

lsqcurvefit stopped because the final change in the sum of squares relative to 
its initial value is less than the value of the function tolerance.

<stopping criteria details>

이 경우에는 이러한 유한 차분 옵션을 사용해도 해가 개선되지는 않습니다.

disp([presnorm,presnorm2])

   20.0637   20.0637

상미분 방정식(ODE) 피팅하기

로렌츠 방정식: 정의 및 수치 해

ODE 해에 원형 경로 피팅하기

원호에 ODE 피팅하기

ODE 피팅의 문제

참고 항목

도움말 항목