最小长度解编程与实践.docVIP

最小长度解编程与实践 班级：地物11101 成员：王云鹤(组长) 吴迪 张林 杨寒 一．小组分工 杨寒和张林涛负责讨论所给题目的思路及程序构图 王云鹤负责编程 吴迪负责整理并完成报告 二．方法原理 对于 （1） 式中：d为M×1维向量；G为 M×N阶矩阵；m为N×1维向量。 1.当MN=r时，方程是超定方程，无常规意义下的解，但有最小方差解。 设e为观测数据d与理论计算值Gm之误差向量，则方差（即目标函数）为: 最小方差解必须满足： 所以 2.当NM=r时，方称是纯欠定方程，也无常规意义上的解，但有最小模型长度。此时的目标函数在（1）式的约束下有极小，即 E=mTm 在（1）式下为极小。 根据极值理论，必须引入拉格朗日算子“λ” 将条件极值问题化为无条件极值问题。因此，目标函数应为： 显然，求上述目标函数极小值问题，可以化为求： 故有 将（5）带入（1），则 再将（6）带入（5），得纯欠定问题的正态方程： 3.当min(M,N)r时，方程是混定方程，也无常规意义下的解，既有超定问题，也有欠定问题的性质。 此时设（1）式是一混定问题，根据混定方程的性质，此时的目标函数应兼有方差项(d-Gm)T(d-Gm)和模型长度项mTm两项内容，即： E=d-GmTd-Gm+ε2mTm =dTd-mTGTd-dTGm+mTGTGm+ε2mTm 求E相对于m或mT的偏导数，并设其为零。简化后得： GTG+ε2Im=GTd 因而 m=(GTG+ε2I)-1GTd 式中ε2为阻尼系数或加权因子。 三．程序框图 超定： 混定： 四．程序代码 1、 超定： program main implicit none integer,parameter::n=3,n1=2 real*8::d(n),g(n,n1),m(n1),gt(n1,n),gtg(n1,n1),gtg1(n1,n1),gtd(n1),gm(n),dgm(n) integer::i,j,l,is(n1),js(n1) d=(/3,1,1/) g(:,1)=(/1,1,0/) g(:,2)=(/1,0,1/) gt=TRANSPOSE (g) gtg=MATMUL(gt,g) gtd=matmul(gt,d) gtg1=gtg call brinv(gtg1,n1,L,IS,JS) m=matmul(gtg1,gtd) write(*,*)m pause end program 2、混定： program main implicit none integer,parameter::n=3,n1=3 real*8::d(n),g(n,n1),m(n1),gt(n1,n),gtg(n1,n1),gtg1(n1,n1),gtd(n1),gm(n),dgm(n),e(n1,n1) real*8::t integer::i,j,l,is(n1),js(n1) d=(/3,1,1/) g(:,1)=(/1,0,0/) g(:,2)=(/1,0,0/) g(:,3)=(/0,1,-1/) e=0 do i=1,n1 e(i,i)=1 enddo gt=TRANSPOSE (g) gtg=MATMUL(gt,g) gtd=matmul(gt,d) t=0 do j=1,20 t=t+0.1 write(*,*)收敛系数,t E=t*t*E GtG1=GtG+E call BRINV(gtg1,N,L,IS,JS) m=matmul(GtG1,gtd) write(*,*)求得结果,m enddo pause endprogram 五．运行结果 超定： 超定方程运算结果： m=4343 与所给答案相符， 混定： 运算结果： m=32320 与所给答案稍有偏差，问题正在查找中。 六．成员签名