AJO_1617_ej4nolineal_triangulo

hace 203 días por etlopez

alfa,beta,gamma = 31.9857, 57.5791, 110.4692 a,b,c = 723.047, 1179.908, 1481.252 sang = 0.01 sdist = 0.001 #tomamos como variables x=alfa, y=beta, z=a F(x,y,z) = (x,y,200-x-y,z,z*sin(pi/200*beta)/sin(pi/200*alfa),z*sin(pi/200*(200-alfa-beta))/sin(pi/200*alfa)) #tomando como vector inicial directamente las observaciones X0 = vector([alfa,beta,a]).column() 
       
X0 
       
[31.9857000000000]
[57.5791000000000]
[723.047000000000]
[31.9857000000000]
[57.5791000000000]
[723.047000000000]
#tomando otro vector inicial X00=vector([200-beta-gamma,200-alfa-gamma,b*sin(alfa*pi/200).n()/sin(beta*pi/200).n()]).column(); X00 
       
[31.9517000000000]
[57.5451000000000]
[722.812864765784]
[31.9517000000000]
[57.5451000000000]
[722.812864765784]
L = vector([alfa, beta, gamma, a,b,c]).column() SigmaL = diagonal_matrix([sang^2,sang^2,sang^2,sdist^2,sdist^2,sdist^2]) #no sabemos cuál es la vdr, tomamos por ejemplo sdist^2 vdr = sdist^2 QL = 1/vdr*SigmaL P = QL^(-1) 
       
show(P) 
       

                                
                            

                                
#podemos ver si los residuos para el vector inicial son altos o no V0 = L - F(*X00.column(0)).column().n() sumaV0 = transpose(V0)*P*V0 sumaV0 = sumaV0[0,0] sumaV0 
       
0.198984566938483
0.198984566938483
# V0 = L - F(*X0.column(0)).column().n() sumaV0 = transpose(V0)*P*V0 sumaV0 = sumaV0[0,0] sumaV0 
       
0.156095761816280
0.156095761816280
#los residuos salen más bajos, indicando que este punto es mejor... 
       
J0 = jacobian(F,(x,y,z))(*X0.column(0)) DeltaL0 = L - F(*X0.column(0)).column() N=transpose(J0)*P*J0 T=transpose(J0)*P*DeltaL0 DeltaX0 = N^(-1) * T DeltaX0 = DeltaX0.n() print "Delta X0:" show(DeltaX0) X1= X0 + DeltaX0 print "Soluciones" show(X0,X1) V1 = L - F(*X1.column(0)).column().n() print "Suma residuos anterior",sumaV0 # la anterior sumaV1 = transpose(V1)*P*V1 sumaV1 = sumaV1[0,0] print "Suma residuos actual",sumaV1 
       
Delta X0:
Soluciones
Suma residuos anterior 0.156095761816280
Suma residuos actual 0.0687062507277703
Delta X0:
Soluciones
Suma residuos anterior 0.156095761816280
Suma residuos actual 0.0687062507277703
#la precisión aún no es la deseada (nos fijamos en los valores de Delta X0), aunque veremos que en realidad X2 será igual que X1 
       
J1 = jacobian(F,(x,y,z))(*X1.column(0)) DeltaL1 = L - F(*X1.column(0)).column() N=transpose(J1)*P*J1 T=transpose(J1)*P*DeltaL1 DeltaX1 = N^(-1) * T DeltaX1 = DeltaX1.n() print "Delta X1:" show(DeltaX1) X2= X1 + DeltaX1 print "Soluciones" show(X0,X1,X2) V2 = L - F(*X2.column(0)).column().n() print "Suma residuos",sumaV0, sumaV1 # la anterior sumaV2 = transpose(V2)*P*V2 sumaV2 = sumaV2[0,0] print "Suma residuos actual",sumaV2 
       
Delta X1:
Soluciones
Suma residuos 0.156095761816280 0.0687062507277703
Suma residuos actual 0.0687062507277703
Delta X1:
Soluciones
Suma residuos 0.156095761816280 0.0687062507277703
Suma residuos actual 0.0687062507277703
#Como podemos ver, el punto ni siquiera varía, nos podemos quedar con X2 o con X1. #Calculemos los errores asociados a estas mediciones, con la varianza de ref a posteriori vdr=(transpose(V2)*P*V2)[0,0]/3 QX=N^(-1) SigmaX=vdr*QX desviaciones=[sqrt(SigmaX[0,0]).n(),sqrt(SigmaX[1,1]).n(),sqrt(SigmaX[2,2]).n()]; print(desviaciones) 
       
[1.23563974188605, 1.23563974188605, 0.0539719384352084]
[1.23563974188605, 1.23563974188605, 0.0539719384352084]
vdr 
       
0.0229020835759234
0.0229020835759234

Los ajustes de las observaciones son:

LmasV = F(*X2.column(0)).n() show(LmasV) 
       

                                
                            

                                
#y la matriz de varianzas asociada a estos ajustes es QLmasV = J1*QX*transpose(J1) QLmasV = QLmasV.n() SigmaLmasV = vdr * QLmasV view(SigmaLmasV) 
       

                                
                            

                                
###PODRÍA TENER SENTIDO HACER UNA PROPAGACIÓN, POR EJEMPLO PARA CALCULAR EL ÁREA DEL TRIÁNGULO 
       
area(x,y,z)=z^2/2*sin((200-x-y)*pi/200)*sin(y*pi/200)/sin(x*pi/200) 
       
area(*X2.column(0)).n() #esta es la estimación del área 
       
420912.963182032
420912.963182032
Jarea=jacobian(area,(x,y,z))(*X2.column(0)).n() 
       
#Ahora, sabemos que la propagación no lineal es, si Y=F(X), QY=Jarea*QX*Jarea^t, SigmaY=Jarea*SigmaX*Jarea^t 
       
SigmaY=Jarea*SigmaX*transpose(Jarea); SigmaY.n() 
       
[3.48543026444556e8]
[3.48543026444556e8]
#la desviación asociada al área es la raíz cuadrada de este número sqrt(SigmaY[0,0]).n() 
       
18669.3070692127
18669.3070692127
 
       
LmasV 
       
(31.9743666666667, 57.5677666666667, 110.457866666667, 722.941575627918,
1180.11810552430, 1481.13605273688)
(31.9743666666667, 57.5677666666667, 110.457866666667, 722.941575627918, 1180.11810552430, 1481.13605273688)
area2(x,y,z)=sqrt((x+y+z)/2*((x+y+z)/2-x)*((x+y+z)/2-y)*((x+y+z)/2-z)) 
       
area2(LmasV[3],LmasV[4],LmasV[5]).n() 
       
420860.310528771
420860.310528771
Jarea2=jacobian(area2,(x,y,z))(LmasV[3],LmasV[4],LmasV[5]).n() 
       
SigmaX=SigmaLmasV[[3..5],[3..5]]; SigmaX 
       
[0.00291297013845392 0.00475508522006814 0.00596798584845895]
[0.00475508522006814 0.00776212400931488 0.00974204332786087]
[0.00596798584845895 0.00974204332786087  0.0122269894281547]
[0.00291297013845392 0.00475508522006814 0.00596798584845895]
[0.00475508522006814 0.00776212400931488 0.00974204332786087]
[0.00596798584845895 0.00974204332786087  0.0122269894281547]
SigmaY=Jarea2*SigmaX*transpose(Jarea2); SigmaY.n() 
       
[3948.80387245003]
[3948.80387245003]
sqrt(SigmaY[0,0]).n() 
       
62.8395088495290
62.8395088495290