Esta página aún no se ha traducido para esta versión. Puede ver la versión más reciente de esta página en inglés.
En este ejemplo se muestra cómo ajustar los parámetros de una ODE a los datos de dos maneras. El primero muestra un ajuste directo de una ruta circular de velocidad constante a una porción de una solución del sistema de Lorenz, una famosa ODE con dependencia sensible de los parámetros iniciales. El segundo muestra cómo modificar los parámetros del sistema de Lorenz para ajustarse a una trayectoria circular de velocidad constante. Puede utilizar el enfoque adecuado para su aplicación como modelo para ajustar una ecuación diferencial a los datos.
El sistema de Lorenz es un sistema de ecuaciones diferenciales ordinarias (ver).El sistema de Lorenz Para constantes reales
Los valores de Lorenz de los parámetros para un sistema sensible son [x(0),y(0),z(0)] = [10,20,10]
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 durante un pequeño intervalo de tiempo, se ve un poco como movimiento circular uniforme. Trace la solución durante 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])
Las ecuaciones de una ruta circular tienen varios parámetros:
Ángulo
Ángulo
RadioR
VelocidadV
Cambio de tiempo 0t0
cambio 3-D en el espaciodelta
En términos de estos parámetros, determine la posición de la ruta circular para veces.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 ruta circular que mejor se ajuste al sistema Lorenz a veces en la solución ODE, utilice.lsqcurvefit
Con el fin de mantener los parámetros en límites razonables, poner límites en los diversos 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.
Trazar los puntos circulares que mejor se ajuste en los tiempos de la solución 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
Ahora modifique los parámetros
Para ello, escriba un archivo de función que tome los parámetros de la ODE Fit y calcule la trayectoria a lo largo de los tiempos.paramfun
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, utilízese para minimizar las diferencias entre la nueva trayectoria de ODE calculada y el arco circular.lsqcurvefit
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 cambió 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 y cambiado por aproximadamente 10%.sigma
beta
Trace 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
Como se describe en, un optimizador puede tener problemas debido al ruido inherente en las soluciones numéricas de ODE.Optimización de una ecuación de simulación o diferencial ordinaria Si sospecha que la solución no es ideal, quizás porque el mensaje de salida o el indicador de salida indica una posible inexactitud, intente cambiar la diferencia finita. En este ejemplo, utilice un tamaño de paso de diferencia finita mayor 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, el uso de estas opciones de diferenciación finita no mejora la solución.
disp([presnorm,presnorm2])
20.0637 20.0637