AJO_1819_Propagacion_ejemplo1

117 days ago by sevillad1819

Ejemplo de propagación no lineal

Problema: dados un acimut y una distancia con sus precisiones (d.t.), calcular las coordenadas del punto resultante y sus precisiones (d.t. de cada coordenada, correlación, y dibujo de elipses de error).

 

Primero escribimos la función que convierte los datos a coordenadas. Como queremos trabajar en grados, pero las funciones trigonométricas están en radianes, dentro metemos el factor de conversión.

f(alfa,d) = ( d*sin(pi*alfa/200) , d*cos(pi*alfa/200) ) f 
       
(alfa, d) |--> (d*sin(1/200*pi*alfa), d*cos(1/200*pi*alfa))
(alfa, d) |--> (d*sin(1/200*pi*alfa), d*cos(1/200*pi*alfa))

A continuación tenemos los datos. Evaluamos la función en esos datos para asegurarnos de que está bien.

ld, sd = 50, 0.01 lalfa, salfa = 70, 0.0020 f(lalfa, ld).n() point(f(lalfa, ld),size=30,aspect_ratio=1,figsize=4) + point((0,0)) 
       

Con las precisiones construimos la matriz de covarianzas inicial. El orden de las variables debe ser el que usamos con la función.

S1 = diagonal_matrix([salfa^2, sd^2]) S1 
       
[ 4.00000000000000e-6    0.000000000000000]
[   0.000000000000000 0.000100000000000000]
[ 4.00000000000000e-6    0.000000000000000]
[   0.000000000000000 0.000100000000000000]

Calculamos el jacobiano y lo guardamos evaluado en los datos.

J = jacobian(f,(alfa,d)) J = J(lalfa,ld) J 
       
[ 1/4*pi*cos(7/20*pi)         sin(7/20*pi)]
[-1/4*pi*sin(7/20*pi)         cos(7/20*pi)]
[ 1/4*pi*cos(7/20*pi)         sin(7/20*pi)]
[-1/4*pi*sin(7/20*pi)         cos(7/20*pi)]

Ahora simplemente hacemos la propagación. La matriz de covarianzas resultante contiene la información estadística de las coordenadas, en el orden en que escribimos la función.

S2 = (J*S1*transpose(J)).n() S2 
       
[0.0000798978121756447 0.0000394527650077175]
[0.0000394527650077175 0.0000225695889246277]
[0.0000798978121756447 0.0000394527650077175]
[0.0000394527650077175 0.0000225695889246277]

A continuación calculamos las desviaciones típicas y la correlación entre las coordenadas.

sx, sy = sqrt(S2[0,0]), sqrt(S2[1,1]) rxy = S2[0,1]/sx/sy sx; sy; rxy 
       
0.00893855761158615
0.00475074614398914
0.929069445421225
0.00893855761158615
0.00475074614398914
0.929069445421225

Pasamos a la elipse de error. En otra hoja tienes el resumen teórico con algún ejemplo de cálculo.

Los semiejes mayor y menor son las raíces cuadradas de los autovalores de la matriz (aunque no sepas nada de autovalores, basta con usar el código siguiente para calcularlos). El semieje mayor es un indicador importante de la precisión.

e = map(sqrt,S2.eigenvalues()) smax = max(e).n() smin = min(e).n() smax; smin 
       
__main__:2: UserWarning: Using generic algorithm for an inexact ring,
which will probably give incorrect results due to numerical precision
issues.
0.0100000000000000
0.00157079632679490
__main__:2: UserWarning: Using generic algorithm for an inexact ring, which will probably give incorrect results due to numerical precision issues.
0.0100000000000000
0.00157079632679490

A continuación se define una función con la que podrás dibujar las elipses que quieras. En la hoja de resumen de elipses de error tienes esta misma función, quizás un poco más actualizada, para usarla en cualquier ejercicio.

def elipse(centro,matriz,prob,escala=1,color="blue"): c = sqrt(r.qchisq(prob,2))._sage_() e = map(sqrt,matriz.eigenvalues()) e = map(abs,e) smax = max(e).n() smin = min(e).n() try: alfa = arctan(matriz[0,1]/(smax^2-matriz[0,0])) except ZeroDivisionError: alfa = 0 return ellipse(centro, escala*c*smax, escala*c*smin, pi/2-alfa, color=color) elipse 
       
<function elipse at 0x7f94821ce2a8>
<function elipse at 0x7f94821ce2a8>

La usamos con nuestro ejemplo para dibujar al 95% y al 99.9%. El centro es la estimación puntual, que se obtiene directamente evaluando la función como hicimos arriba. Es mejor guardar el gráfico en una variable y dibujarlo después con opciones adicionales (tamaño, relación entre los ejes, etc.)

punto = f(lalfa, ld) gelipse1 = elipse(punto, S2, 0.95) gelipse2 = elipse(punto, S2, 0.999, color="red") (gelipse1+gelipse2).show(aspect_ratio=1, figsize=4) 
       

Dibujamos la elipse en referencia al origen de coordenadas.

gorigen = point([0,0]) (gelipse1+gorigen).show(aspect_ratio=1, figsize=7) 
       

Ahora cambia los datos iniciales y ejecuta el código de nuevo para ver cómo cambia la forma de la elipse. Puedes ejecutarlo todo de una vez en el menú "Acción" -> "Evaluar todo".