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, $$\sigma ,\rho ,\beta$$, el sistema es

$$\begin{array}{l}\frac{dx}{dt}=\sigma (y-x)\\ \frac{dy}{dt}=x(\rho -z)-y\\ \frac{dz}{dt}=xy-\beta z.\end{array}$$

Los valores de Lorenz de los parámetros para un sistema sensible son $$\sigma =10,\beta =8/3,\rho =28$$. 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 $$\theta (1)$$ de la trayectoria desde el plano x-y

• Ángulo $$\theta (2)$$ 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 $$\sigma ,\beta ,and\rho$$ 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 salto 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

Temas relacionados

• Discretized Optimal Trajectory, Problem-Based
• Optimizing a Simulation or Ordinary Differential Equation

×

Abrir ejemplo

Tiene una versión modificada de este ejemplo. ¿Desea abrir este ejemplo con sus modificaciones?

No, sobrescribir la versión modificada Sí