Ajustar una ecuación diferencial ordinaria (ODE)

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.

Sistema de Lorenz: definición y solución numérica

El sistema de Lorenz es un sistema de ecuaciones diferenciales ordinarias (ver).El sistema de Lorenz Para constantes reales

<math display="block">

<mrow>

<mi>σ</mi>

<mo>,</mo>

<mspace width="0.5em"></mspace>

<mi>ρ</mi>

<mo>,</mo>

<mspace width="0.5em"></mspace>

<mi>β</mi>

</mrow>

</math>

, el sistema es

<math display="block">

<mrow>

<mtable columnalign="left">

<mtr>

<mtd>

<mrow>

<mfrac>

<mrow>

<mi>d</mi>

<mi>x</mi>

</mrow>

<mrow>

<mi>d</mi>

<mi>t</mi>

</mrow>

</mfrac>

<mo>=</mo>

<mi>σ</mi>

<mo stretchy="false">(</mo>

<mi>y</mi>

<mo>-</mo>

<mi>x</mi>

<mo stretchy="false">)</mo>

</mrow>

</mtd>

</mtr>

<mtr>

<mtd>

<mrow>

<mfrac>

<mrow>

<mi>d</mi>

<mi>y</mi>

</mrow>

<mrow>

<mi>d</mi>

<mi>t</mi>

</mrow>

</mfrac>

<mo>=</mo>

<mi>x</mi>

<mo stretchy="false">(</mo>

<mi>ρ</mi>

<mo>-</mo>

<mi>z</mi>

<mo stretchy="false">)</mo>

<mo>-</mo>

<mi>y</mi>

</mrow>

</mtd>

</mtr>

<mtr>

<mtd>

<mrow>

<mfrac>

<mrow>

<mi>d</mi>

<mi>z</mi>

</mrow>

<mrow>

<mi>d</mi>

<mi>t</mi>

</mrow>

</mfrac>

<mo>=</mo>

<mi>x</mi>

<mi>y</mi>

<mo>-</mo>

<mi>β</mi>

<mi>z</mi>

<mo>.</mo>

</mrow>

</mtd>

</mtr>

</mtable>

</mrow>

</math>

Los valores de Lorenz de los parámetros para un sistema sensible son

<math display="block">

<mrow>

<mi>σ</mi>

<mo>=</mo>

<mn>1</mn>

<mn>0</mn>

<mo>,</mo>

<mspace width="0.5em"></mspace>

<mi>β</mi>

<mo>=</mo>

<mn>8</mn>

<mo>/</mo>

<mn>3</mn>

<mo>,</mo>

<mspace width="0.5em"></mspace>

<mi>ρ</mi>

<mo>=</mo>

<mn>2</mn>

<mn>8</mn>

</mrow>

</math>

. Inicie el sistema desde y visualice la evolución del sistema desde el tiempo 0 hasta el 100.[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])

Ajuste una ruta circular a la solución ODE

Las ecuaciones de una ruta circular tienen varios parámetros:

• Ángulo

  <math display="block">

  <mrow>

  <mi>θ</mi>

  <mo stretchy="false">(</mo>

  <mn>1</mn>

  <mo stretchy="false">)</mo>

  </mrow>

  </math>
  de la ruta del plano x-y

• Ángulo

  <math display="block">

  <mrow>

  <mi>θ</mi>

  <mo stretchy="false">(</mo>

  <mn>2</mn>

  <mo stretchy="false">)</mo>

  </mrow>

  </math>
  del plano desde una inclinación a lo largo del eje x

• 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

Ajuste la oda al arco circular

Ahora modifique los parámetros

<math display="block">

<mrow>

<mi>σ</mi>

<mo>,</mo>

<mspace width="0.5em"></mspace>

<mi>β</mi>

<mo>,</mo>

<mspace width="0.5em"></mspace>

<mrow>

<mstyle mathvariant="normal">

<mrow>

<mi>a</mi>

<mi>n</mi>

<mi>d</mi>

</mrow>

</mstyle>

</mrow>

<mspace width="0.5em"></mspace>

<mi>ρ</mi>

</mrow>

</math>

para ajustarse mejor al arco circular. Para un ajuste aún mejor, permita que el punto inicial [10, 20, 10] cambie también.

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.paramfunt

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.lsqcurvefitsoln

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%.sigmabeta

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

Problemas en el montaje de ODEs

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 

Temas relacionados

• Optimización de una ecuación de simulación o diferencial ordinaria

×

Abrir ejemplo

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

No, sobrescribir la versión modificada Sí