|
|
山中和義さんへのお返事です。
以下プログラム実行しましたが有効桁数が反映されずにフル桁で表示されます。
十進ベーシックの、バージョンが違うからかな？？？
それでも、参考図書を中途半端というか自分の努力分しか読めないのでうれしい限りですが。
山中様どうかお願いですから、プログラムの解説を詳細に書いていただけないでしょうか。
こっちとしても、頑張り不足を痛感しています。
素敵なプログラムをありがとうございました。
> kikiririさんへのお返事です。
>
>
> 参考サイトのURLが変更されているようです。
> 手順および理論は、このサイトのPDFファイルを参照してください。
>
>
> ※「ポケコンプログラムによる測量計算法（山海堂）」の著者のようです。
>
>
> また、前回紹介した本の手計算部分をプログラムしたものを掲載します。参考にしてください。
>
> <PRE>
> !最小２乗法による測量網平均（観測方程式）
>
> !参考文献 最小ニ乗法と測量網平均の基礎 東洋出版 p.170～173（p.165～177）
>
> !既知点
> LET H301=41.7065 !301点の標高（m）
> LET H302=43.9632 !302点
>
>
> !観測値
> LET n=6 !数
>
> DATA 5.1342, 4.6 !路線１ 観測比高（高低差）（m）、路線距離（km）
> DATA -2.0508, 5.2 !路線２
> DATA -2.2287, 9.2 !路線３
> DATA -0.0330, 6.3 !路線４
> DATA -0.8601, 3.9 !路線５
> DATA -0.8258, 8.3 !路線６
>
> DIM Hb(n),P(n,n) !観測値、重量行列
> MAT P=ZER
> FOR i=1 TO n
> READ Hb(i),Wi
> LET P(i,i)=1/Wi
> NEXT i
>
> MAT PRINT Hb; !debug
> MAT PRINT P; !debug
>
>
> !新点（未知点）
> LET m=3 !数
> LET pnt$="１２３"
>
>
> !------------------------------ 観測方程式の一般式 V=A*X+Lを組み立てる p.170～171
>
> DIM A(n,m) !係数行列 ※各路線の観測比高と新点標高との関係より
> DATA -1, 1, 0
> DATA 0,-1, 0
> DATA 0, 0, 1
> DATA 0, 0,-1
> DATA 1, 0,-1
> DATA 1, 0, 0
> MAT READ A
>
> DIM L(n,1) !定数項
> LET L(1,1)=-Hb(1)
> LET L(2,1)=-(Hb(2)-H302)
> LET L(3,1)=-(H302+Hb(3))
> LET L(4,1)=-(Hb(4)-H301)
> LET L(5,1)=-Hb(5)
> LET L(6,1)=-(H301+Hb(6))
>
> MAT PRINT A; !debug
> MAT PRINT L; !debug
>
>
> !------------------------------ 方程式の解（X~）を求める p.172中
>
> DIM TT(50,50),MM(50,50) !作業用
>
> MAT TT=TRN(A) !TRN(A)*P*Aの計算
> MAT TT=TT*P
> MAT MM=TT !save it
> MAT TT=TT*A
> MAT PRINT TT; !debug
>
> DIM Q(m,m) !重み係数行列 Q=INV(TRN(A)*P*A)の計算
> MAT Q=INV(TT)
> MAT PRINT Q; !debug
>
>
> MAT TT=MM*L !TRN(A)*P*Lの計算
> MAT PRINT TT; !debug
>
>
> DIM Xb(m,1) !最確値 X~
> MAT TT=Q*TT
> MAT Xb=(-1)*TT
> MAT PRINT Xb; !debug
>
>
> !------------------------------ 網平均の最終結果 p.172中～173
>
> DIM Vb(n,1) !残差ベクトル V~=L-A*INV(TRN(A)*P*A)*TRN(A)*P*Lの計算
> MAT TT=A*Xb !Xb=-INV(TRN(A)*P*A)*TRN(A)*P*Lより
> MAT Vb=L+TT
> MAT PRINT Vb; !debug
>
> MAT TT=TRN(Vb) !単位重みの標準偏差σ0の推定値 m0=SQR(TRN(V~)*P*V~/(m-n))の計算
> MAT TT=TT*P
> MAT TT=TT*Vb
> MAT PRINT TT; !debug
> LET m0=SQR(TT(1,1)/(n-m))
> PRINT "m0="; m0 !debug
>
> FOR i=1 TO m
> PRINT pnt$(i:i); "点の標高"; Xb(i,1);"m"; !有効桁数 ##.####
> PRINT " ±"; m0*SQR(Q(i,i)); "mm" !有効桁数 #.#
> NEXT i
>
>
> END
> </PRE>
|
|