AJO_1617_Propagación_ejemplo2

hace 187 días por sevillad1617

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 = 30,0.02 lb,sb = 40,0.02 lc,sc = 50,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() 
       
(40.9665529398267, 59.0334470601733, 100.000000000000)
(40.9665529398267, 59.0334470601733, 100.000000000000)

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.0371209926798273, 0.0543514065930112, 0.0750263596797588]
[0.0371209926798273, 0.0543514065930112, 0.0750263596797588]

Así podemos calcular el intervalo al 95% para el ángulo recto.

print f1(la,lb,lc)[2],"+-",2*sqrt(S21[2,2]) 
       
100 +- 0.150052719359518
100 +- 0.150052719359518

Probamos a cambiar la función por otra que se nos podría haber ocurrido; por ejemplo, en el tercer ángulo, en vez del teorema del coseno usaremos que es el complementario a 200 de los otros dos. Veamos que la nueva 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() 
       
(40.9665529398267, 59.0334470601733, 100.000000000000)
(40.9665529398267, 59.0334470601733, 100.000000000000)
(40.9665529398267, 59.0334470601733, 100.000000000000)
(40.9665529398267, 59.0334470601733, 100.000000000000)

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() 
       
(40.9665529398267, 59.0334470601733, 100.000000000000, 50.0000000000000)
(40.9665529398267, 59.0334470601733, 100.000000000000, 50.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, todo en la última fila y la última columna, que son las nuevas.

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

                                
                            

                                

La correlación entre el ángulo $C$ (tercera componente de la función) y el lado $c$ (cuarta componente) se calcula fácilmente. La primera línea del bloque siguiente calcula la lista de desviaciones típicas.

s = map(sqrt,S23.diagonal()) S23[2,3]/s[2]/s[3] 
       
0.707106781186548
0.707106781186548

La correlación es positiva, lo que tiene sentido: cuanto más grande un lado, más grande su ángulo opuesto. Calculamos las otras dos correlaciones (es un buen ejercicio de geometría el preguntarse por qué salen esos signos).

S23[0,3]/s[0]/s[3]; S23[1,3]/s[1]/s[3] 
       
-0.514495755427527
-0.624695047554424
-0.514495755427527
-0.624695047554424