AJO_1819_Propagación_ejemplo2

118 days ago by sevillad1819

Ejemplo de propagación no lineal

Problema: dados los lados de un triángulo con sus precisiones (d.t.), calcular los tres ángulos y sus desviaciones típicas.

 

Primero escribimos la función que calcula cada ángulo. En los tres casos usamos el teorema del coseno, debajo lo haremos de otra manera. Como queremos trabajar en grados, pero las funciones trigonométricas están en radianes, dentro metemos el factor de conversión.

la,sa = 15,0.02 lb,sb = 20,0.02 lc,sc = 30,0.02 f1(a,b,c) = ( 200/pi*arccos((a^2-b^2-c^2)/(-2*b*c)), 200/pi*arccos((b^2-a^2-c^2)/(-2*a*c)), 200/pi*arccos((c^2-b^2-a^2)/(-2*b*a)) ) f1(la,lb,lc).n() 
       
(29.3159219437866, 40.3733972384599, 130.310680817753)
(29.3159219437866, 40.3733972384599, 130.310680817753)

Hacemos la propagación y obtenemos la matriz 3x3 que contiene la información sobre la precisión de los ángulos y sus relaciones.

S1 = diagonal_matrix([sa^2,sb^2,sc^2]) J1 = jacobian(f1,(a,b,c)) J1 = J1(la,lb,lc) S21 = (J1*S1*transpose(J1)).n() view(S21) 
       

                                
                            

                                

Con el comando map aplicamos una función (en este caso la raíz cuadrada) a una lista (en este caso la diagonal).

map(sqrt,S21.diagonal()) 
       
[0.0976616594378941, 0.135488275013595, 0.224298330207189]
[0.0976616594378941, 0.135488275013595, 0.224298330207189]

En este caso las correlaciones no nos interesan.

 

Ahora cambiamos la fórmula del tercer ángulo. Veamos que la función produce el mismo resultado que antes.

f2(a,b,c) = ( 200/pi*arccos((a^2-b^2-c^2)/(-2*b*c)), 200/pi*arccos((b^2-a^2-c^2)/(-2*a*c)), 200-200/pi*arccos((a^2-b^2-c^2)/(-2*b*c))-200/pi*arccos((b^2-a^2-c^2)/(-2*a*c)) ) f1(la,lb,lc).n(); f2(la,lb,lc).n() 
       
(29.3159219437866, 40.3733972384599, 130.310680817753)
(29.3159219437866, 40.3733972384599, 130.310680817753)
(29.3159219437866, 40.3733972384599, 130.310680817753)
(29.3159219437866, 40.3733972384599, 130.310680817753)

$A$ y $c$ son la primera y cuarta componente de la función.

Hacemos la propagación con la nueva función. La nueva matriz es la misma que la antigua.

J2 = jacobian(f2,(a,b,c)) J2 = J2(la,lb,lc) S22 = (J2*S1*transpose(J2)).n() view(S21); print; view(S22) 
       


Para terminar, supongamos que queremos ver cómo afecta la variabilidad del lado $c$ a los ángulos. Para eso vamos a calcular las correlaciones. Incluimos el lado en la función como una componente más:

f3(a,b,c) = ( 200/pi*arccos((a^2-b^2-c^2)/(-2*b*c)), 200/pi*arccos((b^2-a^2-c^2)/(-2*a*c)), 200/pi*arccos((c^2-b^2-a^2)/(-2*b*a)), c ) f3(la,lb,lc).n() 
       
(29.3159219437866, 40.3733972384599, 130.310680817753, 30.0000000000000)
(29.3159219437866, 40.3733972384599, 130.310680817753, 30.0000000000000)

Hacemos la propagación. Fíjate en que la matriz 3x3 de arriba es una submatriz de esta. Ahora simplemente tenemos información adicional: la varianza de $c$ (aunque ya la teníamos) y las covarianzas con los tres ángulos.

J3 = jacobian(f3,(a,b,c)) J3 = J3(la,lb,lc) S23 = (J3*S1*transpose(J3)).n() view(S23) 
       

                                
                            

                                
s = map(sqrt,S23.diagonal()) S23[0,3]/s[0]/s[3] 
       
-0.590822551825122
-0.590822551825122

La correlación es negativa, lo que tiene sentido: si aumenta el lado $c$ disminuye el ángulo adyacente $A$. Calculamos las otras dos.

S23[1,3]/s[1]/s[3]; S23[2,3]/s[2]/s[3] 
       
-0.631465906670717
0.638689272171371
-0.631465906670717
0.638689272171371