Ajustar una ecuación diferencial ordinaria (ODE)
Este ejemplo muestra cómo ajustar parámetros de una ODE a datos de dos formas. La primera muestra un ajuste sencillo de una trayectoria circular con velocidad constante a una porción de una solución del sistema de Lorenz, una ODE famosa con dependencia sensible de los parámetros iniciales. La segunda muestra cómo modificar los parámetros del sistema de Lorenz para ajustar una trayectoria circular con velocidad constante. Puede utilizar el enfoque adecuado a la aplicación como un modelo para ajustar una ecuación diferencial a datos.
Sistema de Lorenz: definición y solución numérica
El sistema de Lorenz es un sistema de ecuaciones diferenciales ordinarias (consulte Sistema de Lorenz). En constantes reales, , el sistema es
Los valores de Lorenz de los parámetros para un sistema sensible son . Empiece el sistema desde [x(0),y(0),z(0)] = [10,20,10]
y observe la evolución del sistema del tiempo 0 al 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])
La evolución es bastante complicada, pero, en un intervalo de tiempo pequeño, en cierto modo, parece un movimiento circular uniforme. Represente la solución sobre el intervalo de tiempo [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])
Ajustar una trayectoria circular a la solución de la ODE
Las ecuaciones de una trayectoria circular tienen varios parámetros:
Ángulo de la trayectoria desde el plano x-y
Ángulo del plano desde una inclinación a lo largo del eje x
Radio R
Velocidad V
Desplazamiento t0 desde el tiempo 0
Desplazamiento 3D en delta de espacio
En términos de estos parámetros, determine la posición de la trayectoria circular para los tiempos 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);
Para encontrar la trayectoria circular que mejor se ajusta al sistema de Lorenz en los tiempos que indica la solución de la ODE, utilice lsqcurvefit
. Para mantener los parámetros en límites razonables, ponga límites en los distintos parámetros.
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.
Represente los puntos circulares que mejor se ajustan a los tiempos de la solución de la ODE junto con la solución del sistema de Lorenz.
soln = a + residual; hold on plot3(soln(:,1),soln(:,2),soln(:,3),'r') legend('ODE Solution','Circular Arc') hold off
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
Ajustar la ODE al arco circular
Ahora, modifique los parámetros para que se adapten lo mejor posible al arco circular. Para que se adapten incluso mejor, permita que el punto inicial [10,20,10] también cambie.
Para ello, escriba un archivo de función paramfun
que tome los parámetros del ajuste de la ODE y calcule la trayectoria sobre los tiempos t
.
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);
Para encontrar los mejores parámetros, utilice lsqcurvefit
para minimizar las diferencias entre la nueva trayectoria calculada de la ODE y el arco circular 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.
Determine cuánto ha cambiado esta optimización los parámetros.
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
Los parámetros sigma
y beta
han cambiado en un 10%.
Represente la solución modificada.
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
Problemas al ajustar las ODE
Como se describe en Optimizing a Simulation or Ordinary Differential Equation, un optimizador puede tener problemas debido al ruido inherente a las soluciones numéricas de la ODE. Si sospecha que su solución no es ideal, tal vez porque el mensaje de salida o el indicador de salida indican una posible imprecisión, pruebe a cambiar la diferenciación finita. En este ejemplo, utilice un tamaño de paso de diferencia finita más grande y diferencias finitas centrales.
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.
En este caso, utilizar estas opciones de diferenciación finita no mejora la solución.
disp([presnorm,presnorm2])
20.0637 20.0637