[Editar: Aparentemente, la pregunta implica que hay más de dos puntos disponibles y leí mal la pregunta, por lo que mi enfoque supone solo dos puntos. La pregunta de dos puntos subdeterminada es más divertida, así que dejo esto de todos modos.]
Estas son las ecuaciones clásicas de Lotka-Volterra / Predator-Prey . Asumo la forma:
[matemáticas] \ dot {x} = a_1 x – a_2 xy [/ matemáticas], (1.a)
[matemáticas] \ dot {y} = -b_1 y + b_2 xy [/ matemáticas] (1.b)
La pregunta indica que solo poseemos dos puntos, [matemática] (x_0, y_0) [/ matemática] en [matemática] t_0 [/ matemática] y [matemática] (x_1, y_1) [/ matemática] en [matemática] t_1 [ / math], en el espacio de fase del sistema. Deje que [math] S \ equiv [a_1, a_2, b_1, b_2] ^ T [/ math] sea el conjunto de soluciones deseadas.
El sistema de ecuaciones se puede describir como:
[matemáticas]
\izquierda[
\ begin {matrix}
a_1 x_ {0} & – a_2 x_0 y_0 & + b_1 \ cdot 0 & + b_2 \ cdot 0 \\
a_1 \ cdot 0 & + a_2 \ cdot 0 & – b_1 y_0 & + b_2 x_0 y_0 \\
a_1 x_ {1} & – a_2 x_1 y_1 & + b_1 \ cdot 0 & + b_2 \ cdot 0 \\
a_1 \ cdot 0 & + a_2 \ cdot 0 & – b_1 y_1 & + b_2 x_1 y_1
\ end {matriz}
\ right] =
\izquierda[
\ begin {matrix}
\ dot {x} (t_0) \\
\ dot {y} (t_0) \\
\ dot {x} (t_1) \\
\ dot {y} (t_1)
\ end {matriz}
\derecho]
[/ matemáticas] (2)
Si se nos permitiera obtener información sobre las tasas del sistema ([matemática] \ dot {x} (t) [/ matemática], [matemática] \ dot {y} (t) [/ matemática]) en las dos mediciones del sistema , [matemática] AS = B [/ matemática] sería trivial de resolver por inversión donde,
[matemáticas]
A = \ left [
\ begin {matrix}
x_ {0} y -x_0 y_0 y 0 y 0 \\
0 y 0 y -y_0 y x_0 y_0 \\
x_ {1} y -x_1 y_1 y 0 y 0 \\
0 y 0 y -y_1 y x_1 y_1
\ end {matriz}
\derecho],
[/matemáticas]
[matemáticas]
B = \ left [
\ begin {matrix}
\ dot {x} (t_0) \\
\ dot {y} (t_0) \\
\ dot {x} (t_1) \\
\ dot {y} (t_1)
\ end {matriz}
\derecho],
[/matemáticas]
y [matemáticas] S = A \ barra diagonal inversa B [/ matemáticas]. Pero no tenemos valores de tasa. Entonces no podemos usar esto.
- Análisis numérico: ¿Cuáles son las situaciones en el mundo físico donde surgen problemas de perturbación singulares?
- ¿Qué son las ecuaciones de Lagrange y Hamilton?
- ¿Cuáles son los gastos generales computacionales involucrados al resolver ecuaciones diferenciales parciales en entornos HPC?
- ¿Cuál es la relación y / o diferencia entre un campo de pendiente y un campo vectorial?
- ¿Cuáles son las ventajas de usar un método explícito (por ejemplo, Forward Euler) sobre métodos implícitos, cuando el sistema de ecuaciones diferenciales le gusta esto, [matemática] C dx (t) / dt = G x (t) + u (t) [/ matemática] (donde [matemática] C [/ matemática] y [matemática] G [/ matemática] son matrices dispersas más grandes, [matemática] u (t) [/ matemática] se ingresa en este sistema, [matemática] x (t) [/ matemáticas] es el estado del sistema)?
Por lo tanto, nuestro único recurso es utilizar la información en el flujo (dinámica) del sistema para obtener algunos valores de [math] S [/ math].
Cuando usamos la dinámica del sistema, el problema se convierte en un ejercicio de ajuste de curvas con datos insuficientes ( sistema subdeterminado ) que, en el mejor de los casos, puede ser consistente y nos da un conjunto infinito de resultados ( infinitas maneras de vincular dos puntos en el espacio de fase usando rigidez, no lineal ecuaciones ). Así que reconciliemos con ninguna solución [math] S [/ math] sea única . El sistema Lotka-Volterra también es extremadamente sensible al acondicionamiento numérico, por lo que, a menos que preste especial atención a la precisión de coma flotante, ni siquiera podrá recrear una simulación numérica. La pregunta pide una solución ‘aproximada’, por lo que dejaré de azotar este punto.
El enfoque que presento aquí es tratar este problema como un problema de optimización local restringido estándar donde las trayectorias de fase se obtienen de una solución numérica a las ecuaciones de estado de Lotka-Volterra. Necesitamos un conjunto de conjeturas iniciales para que [math] S [/ math] comience la búsqueda de vecindario. La métrica de la función de costo exacta que se minimiza será un error al cuadrado,
[math] \ text {Error} = \ Sigma_ {i \ in \ {0,1 \}} (x_i ^ p – x_i) ^ 2 + (y_i ^ p – y_i) ^ 2 [/ math], donde [math ] q ^ p [/ math] son los valores predichos / ajustados. Elegí esta métrica para hacer un punto sobre la declaración del problema mal condicionado (ver Fig. 1 a continuación). Puede elegir una métrica diferente según sus creencias.
Conjeturas iniciales
Las conjeturas iniciales localizan el vecindario de la [matemática] S [/ matemática] final. Los nuevos supuestos para las conjeturas iniciales aislarán las nuevas topologías del sistema, cada sistema tendrá una dinámica diferente. Esbozaré algunos enfoques para establecer vecindarios. El enfoque que elija dependerá del proceso exacto de interacción con la población que esté tratando de modelar y sus creencias con respecto a la información inicial que posee. Mezcla y combina y mejora.
La ecuación 2 muestra el sistema de ecuaciones que estaremos resolviendo para obtener [matemáticas] S [/ matemáticas] basadas en creencias. Tenemos una [matemática] A [/ matemática] rígidamente definida. Pero se nos permiten creencias sobre el carácter dinámico del sistema, sobre las tasas [matemáticas] B [/ matemáticas] en el sistema y también se nos permiten nociones preconcebidas sobre [matemáticas] S [/ matemáticas].
Creencias sobre características dinámicas.
- Puntos de equilibrio : este sistema de ecuaciones tiene un nodo de silla de montar en [matemáticas] (x, y) = (0,0) [/ matemáticas] y tiene una constante de movimiento en [matemáticas] (x, y) = (b_1 / b_2 , a_1 / a_2) [/ math] (ciclos inestables centrados en el punto crítico). Podemos elegir cambiar la ubicación del nodo de silla de montar, pero eso no es interesante. También podemos elegir el centro (locus) de un ciclo, más interesante. Llamemos al par ordenado, [matemática] (1 / v, 1 / u) [/ matemática] de modo que [matemática] b_2 = v b_1 [/ matemática] y [matemática] a_2 = u a_1 [/ matemática] (Intente asegúrese de que [matemáticas] (1 / v, 1 / u) [/ matemáticas] no sean irracionales o trascendentales para evitar errores de truncamiento de coma flotante, lo mejor es usar fracciones decimales finitas).
[matemáticas] a_1 = (\ dot {x} _0 / x_0) / (1-uy_0) [/ matemáticas],
[matemáticas] a_2 = u (\ dot {x} _0 / x_0) / (1-uy_0) [/ matemáticas],
[matemáticas] b_1 = (\ dot {y} _0 / y_0) / (1-vx_0) [/ matemáticas],
[matemática] b_2 = v (\ dot {y} _0 / y_0) / (1-vx_0) [/ matemática].
Ahora necesita suposiciones sobre las tasas iniciales [matemáticas] (\ dot {x} _0, \ dot {y} _0) [/ matemáticas]. También puede optar por asumir las tasas finales que se analizan a continuación. Tu decides. - Período aproximado del ciclo : está dado por [math] T = 2 \ pi / \ sqrt {a_1 b_1} [/ math]. Esto deja tres incógnitas en [matemáticas] S [/ matemáticas], y aún necesita creencias sobre las tasas.
- Valores propios : esto equivale a seleccionar [matemática] a_1 [/ matemática] y [matemática] -b_1 [/ matemática] para el valor propio en el nodo de silla de montar. El otro conjunto de valores propios para ciclos viene dado por [math] \ pm \ sqrt {a_1 b_1} [/ math] que es equivalente a establecer el período del ciclo.
- Isoclinas : las dos isoclinas en este sistema son [matemática] (0, a_1 / a_2) [/ matemática] para flujo vertical y [matemática] (b_1 / b_2,0) [/ matemática] para flujo horizontal.
Creencias sobre las tasas, [matemáticas] B [/ matemáticas]
- Elija las tasas iniciales para que sean cero, [matemática] \ dot {x} _0 = 0 [/ matemática] y [matemática] \ dot {y} _0 = 0 [/ matemática] y las tasas finales se aproximarán por una línea, [matemática ] \ dot {x} _1 = (x_1-x_0) / (t_1-t_0) [/ math] y [math] \ dot {y} _1 = (y_1-y_0) / (t_1-t_0) [/ math].
- Elija capturar un ciclo completo. Puede elegir [matemáticas] t_1 – t_0 = T [/ matemáticas] por diversión (use valores aproximados que eviten singularidades, esto tomará una cantidad de tiempo decente para resolver). Necesita [matemática] x_0 \ aprox. X_1 [/ matemática], [matemática] y_0 \ aprox. Y_1 [/ matemática], [matemática] \ dot {x} _0 \ aprox \ dot {x} _1 [/ matemática], [matemática ] \ dot {y} _0 \ aprox \ dot {y} _1 [/ math].
Creencias sobre [matemáticas] S [/ matemáticas]
- Esquema de Las Vegas: fuerza bruta, enfoque de suerte ciega. Adivina [matemáticas] S [/ matemáticas]. Si obtienes una solución, genial. De lo contrario, adivina de nuevo. Esto refleja una creencia totalmente preconcebida en [matemáticas] S [/ matemáticas].
- Preconcepción acerca de un solo elemento en [matemáticas] S [/ matemáticas] – Supongamos que se conoce [matemáticas] a_1 [/ matemáticas], tamaños de [matemáticas] A [/ matemáticas], [matemáticas] B [/ matemáticas] y [matemáticas] S [/ matemáticas] ahora son 3 × 3, 3 × 1, 3 × 1 resp. Por ejemplo, las siguientes matrices suponen que [math] a_1 [/ math] es conocida y las tasas iniciales son 0, mientras que las tasas finales son pendientes lineales entre [math] q_0 [/ math] y [math] q_1 [/ math].
[matemáticas] A = \ left [
\ begin {matrix}
x_0 y_0 & -y_0 y x_0 y_0 \\
-x_1 y_1 y 0 y 0 \\
0 y -y_1 y -x_1 y_1
\ end {matriz}
\ right], [/ math]
[matemáticas] B = \ left [
\ begin {matrix}
a_1 x_0 \\
(x_1-x_0) / (t_1 -t_0) -a_1 x_1 \\
(y_1 -y_0) / (t_1 – t_0)
\ end {matriz}
\ right]. [/ math] - Preconcepción sobre elementos múltiples en [math] S [/ math] – Hágalo usted mismo.
En el script que sigue, primero localizo el vecindario de la solución usando conjeturas iniciales. Luego lo resuelvo restringiendo las tasas del sistema. Las limitaciones se describieron anteriormente. No asumo nada sobre la proximidad de [math] (x_0, y_0) [/ math] y [math] (x_1, y_1) [/ math] a los puntos fijos de este sistema (debe verificar el condicionamiento de la solución )
Aquí hay una figura. 
Figura 1. La figura muestra por qué la pregunta no es tan buena. Muchas formas de ir entre dos puntos en el espacio de fase. Muestro 5 trayectorias (líneas discontinuas) que conectan (x0, y0) y (x1, y1) obtenidas de 5 conjuntos diferentes de supuestos. Hay casos en los que puede conducir el sistema a través de múltiples ciclos antes de alcanzar (x1, y1) como se muestra en la Fig. 2. Utilizo log (x) para comprimir el eje x para la representación visual.
Figura 2. Es posible obtener el mejor ajuste en un ciclo posterior. Probablemente significa condicionamiento de mierda, no otra cosa.
Aquí está el guión.
función LotkaVolterraDataFitting3
cierra todo; clc;
% eliminar (gcp); parpool ('SSHClusterA_N32')
%% Variables de usuario
t0 = 0; t1 = 0,6038;
x0 = 0.118; y0 = 0.2035; x1 = 0.2212; y1 = 1.0235;
T = [t0 t1]; % Del vector que especifica la unidad de tiempo de mediciones, tf = 1 normaliza el período
xini = [x0 x1]; % En t = t0, x = x0; En t = t1, x = x1;
yini = [y0 y1]; % En t = t0, y = y0; En t = t1, y = y1;
%% The 0.5 * Juego de adivinanzas evaluado
% Te muestro un par de formas de hacerlo. Comente / descomente según sea necesario.
% Restringir gradientes de Jacobian, no se muestra en discusión
A = [x0 -x0 * y0 0 0; -y0 x0 * y0 0 0; 0 0 x1 0; 0 0 -y1 x1 * y1];
B = [0; 0; ((x1-x0) / (t1-t0)) + (- (x0 / y0) * x1 * y1); (y1-y0) / (t1-t0)];
S_LousyGuess = A \ B; % Llamándolo como es
% De tasas calculadas de vuelta para mostrar un ciclo completo
% A = [x0 -x0 * y0 0 0; 0 0 -y0 x0 * y0; x1 -x1 * y1 0 0; 0 0 -y1 x1 * y1];
% B = [0.9 -0.5 -2.2 5.4] ';
% S_LousyGuess = A \ B;
% De valor preconcebido para a1
% a1 = -10; % Comencemos con esta suposición, simplemente forzamos un vecindario en a1 manualmente
% A = [x0 * y0 -y0 x0 * y0; -x1 * y1 0 0; 0 -y1 -x1 * y1];
% S_LousyGuess = A \ [a1 * x0; ((x1-x0) / (t1-t0)) - (a1 * x1); (y1-y0) / (t1-t0)];
% S_LousyGuess = [a1 S_LousyGuess ']';
% De valor preconcebido para a1
% a1 = -1;
% A = [x0 * y0 -y0 x0 * y0; -x1 * y1 0 0; 0 -y1 -x1 * y1];
% S_LousyGuess = A \ [a1 * x0; ((x1-x0) / (t1-t0)) - (a1 * x1); (y1-y0) / (t1-t0)];
% S_LousyGuess = [a1 S_LousyGuess ']';
% Forzar diferentes tasas cambiando los límites de tiempo
% t0 = 0; t1 = 0.8;
% T = [t0 t1];
% A = [x0 -x0 * y0 0 0; 0 0 -y0 x0 * y0; x1 -x1 * y1 0 0; 0 0 -y1 x1 * y1];
% B = [0.9 -0.5 -2.2 5.4] ';
% S_LousyGuess = A \ B;
%% Neighborhood búsqueda de error de mínimos cuadrados
disp (['Suposición pésima para S = [' num2str (S_LousyGuess ')'] '])
% Empaquetar todos los datos conocidos en una matriz Sc para su optimización. Tenga en cuenta que esto incluye valores en t = 0
Sc = [xini (1) yini (1) S_LousyGuess '];
[Sc, fval] = fminsearch (@ minimizarLSQErr, Sc, optimset ...
('MaxFunEvals', 10000, 'MaxIter', 10000), T, xini, yini); % De error de búsqueda mínima local
a1 = Sc (3); a2 = Sc (4); b1 = Sc (5); b2 = Sc (6);
disp (['a1 =' num2str (a1) ', a2 =' num2str (a2) ', b1 =' num2str (b1) ', b2 =' num2str (b2)])
%% Recreando el espacio de fase
[~, xm] = ode23 (@ statepace, T ', [xini (1) yini (1) a1 a2 b1 b2]');
plot (xini, yini, 'co', 'MarkerSize', 11, 'LineWidth', 3), grid on, hold on
plot (xm (:, 1), xm (:, 2), 'k -', 'LineWidth', 3), title ('Comprueba esta mierda', 'FontSize', 15, 'FontWeight', 'Bold '),
xlabel ('x', 'FontSize', 15, 'FontWeight', 'Bold'), ylabel ('y', 'FontSize', 15, 'FontWeight', 'Bold'), leyenda ('Datos dados', ' Trayectoria estimada ')
fin
función LSQErr = minimizarLSQErr (coeff, tval, xval, yval)
n = longitud (tval);
[~, y] = ode23 (@ SearchStateSpace, tval, [coeff (1), coeff (2)], [], coeff (3), coeff (4), coeff (5), coeff (6));
epsilonx = y ([1 extremo], 1) -xval (1: n) ';
epsilony = y ([1 final], 2) -yval (1: n) ';
LSQErr = epsilonx '* epsilonx + epsilony' * epsilony; % De error métrica
fin
función dx = espacio de estado (~, y)
dx = [0; 0];
x = [y (1); y (2)];
a1 = y (3); a2 = y (4); b1 = y (5); b2 = y (6); % <-Ser perezoso.
dx (1) = a1 * x (1) - a2 * x (1) * x (2);
dx (2) = - b1 * x (2) + b2 * x (1) * x (2);
dx (3) = 0; dx (4) = 0; dx (5) = 0; dx (6) = 0; % <-Ser perezoso, con cuidado.
fin
función dy = SearchStateSpace (~, y, a1, a2, b1, b2)
dy = [0; 0];
dy (1) = a1 * y (1) - a2 * y (1) * y (2);
dy (2) = -b1 * y (2) + b2 * y (1) * y (2);
fin
(Título de la figura: respuesta de Jacob VanWagoner a ¿Cuál es la cosa más divertida que has visto en una tarea de clase de física?)
YMMV. No use valores iniciales como [math] \ pi [/ math]. Juegue con algunos de los valores iniciales pero no use factores de 100, sea razonable; no asuma un conjunto de espacios de fase como [x0=0.01; x1=0.0102; y0=0.001; y1=10.2;] [x0=0.01; x1=0.0102; y0=0.001; y1=10.2;] [x0=0.01; x1=0.0102; y0=0.001; y1=10.2;] sin modificar las condiciones de contorno apropiadas. Las ecuaciones de Lotka-Volterra se pueden hacer muy rígidas, pero los comportamientos rígidos NO se describen con estas ecuaciones. Así que no hagas rígida tu trayectoria, el código se torcerá. Si lo desea, verifique los valores propios del sistema linealizado jacobiano y su proximidad a las condiciones iniciales. etc.
Este código puede extenderse (con modificación de notación) para su uso en el ajuste de mínimos cuadrados real de datos medidos a los sistemas de tipo de presa Lotka-Volterra / Predator específicos o generalizados (ecuaciones ópticas de Bloch, sistemas de interacción superior asimétrica …) y puede usarse para obtener ajustes reales y [matemáticas] R ^ 2 [/ matemáticas] / pruebas de bondad de ajuste. Y no necesita suposiciones 0.5Assed para un sistema sobredeterminado. Aún necesita redefinir su métrica de error para tener en cuenta todas las mediciones (no el primer y último punto como el que uso aquí).