matlab 常见经典平差 程序
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
希望本人编写的经典平差可以对matlab初学者以及测绘专业的学生有用
By cowboy
function void()
sprintf('请选择平差类型')
sprintf('1: 参数平差\n2: 条件平差\n3：具有条件的参数平差\n4: 具有参数的条件平差\n') a=input('请输入对应号码\n');
switch (a)
case 1
% 参数平差V=AX-L
sprintf('1:参数平差V=AX-L')
n=input('请输入n=');
t=input('请输入t=');
A=input('请输入系数矩阵A=');
L=input('请输入常数项矩阵L=');
P=diag(input('请输入权阵P='))
%P=input('请输入权阵P=');
N=A'*P*A;
U=A'*P*L;
sprintf('开始计算')
X=inv(N)*(U);
V=A*X-L;
Qx=inv(N);
u=sqrt((V'*P*V)/(n-t));
fid=fopen('data.txt','wt'); %写入文件路径，文件输出
fprintf(fid,'=====================参数平差: V=AX-L =====================\n');
fprintf(fid,'输出n,t\n') ;
fprintf(fid,'%12.5f %12.5f\n',n ,t) ;%n,t
fprintf(fid,'输出矩阵A\n')
[m,n]=size(A); %A
for i=1:1:m
for j=1:1:n
if j==n
fprintf(fid,'%12.5f\n',A(i,j));
else
fprintf(fid,'%12.5f\t',A(i,j));
end
end
end
fprintf(fid,'\n');
fprintf(fid,'输出矩阵L\n') ;
[m,n]=size(L); %L
for i=1:1:m
for j=1:1:n
if j==n
fprintf(fid,'%12.5f\n',L(i,j));
else
fprintf(fid,'%12.5f\t',L(i,j));
end
end
end
fprintf(fid,'\n');
fprintf(fid,'输出矩阵P\n') ;
[m,n]=size(P); %P for i=1:1:m
for j=1:1:n
if j==n
fprintf(fid,'%12.5f\n',P(i,j));
else
fprintf(fid,'%12.5f\t',P(i,j));
end
end
end
fprintf(fid,'\n');
fprintf(fid,'输出矩阵X\n');
[m,n]=size(X); %X for i=1:1:m
for j=1:1:n
if j==n
fprintf(fid,'%12.5f\n',X(i,j));
else
fprintf(fid,'%12.5f\t',X(i,j));
end
end
end
fprintf(fid,'\n');
fprintf(fid,'输出矩阵V\n');
[m,n]=size(V); %V for i=1:1:m
for j=1:1:n
if j==n
fprintf(fid,'%12.5f\n',V(i,j));
else
fprintf(fid,'%12.5f\t',V(i,j));
end
end
end
fprintf(fid,'\n');
fprintf(fid,'输出矩阵Qx\n');
[m,n]=size(Qx); %Qx for i=1:1:m
for j=1:1:n
if j==n
fprintf(fid,'%12.5f\n',Qx(i,j));
else
fprintf(fid,'%12.5f\t',Qx(i,j));
end
end
end
fprintf(fid,'\n');
fprintf(fid,'输出单位权中误差u\n');
[m,n]=size(u); %u for i=1:1:m
for j=1:1:n
if j==n
fprintf(fid,'%12.5f\n',u(i,j));
else
fprintf(fid,'%12.5f\t',u(i,j));
end
end
end
fprintf(fid,'\n');
fclose(fid);
case 2
% 条件平差BV+W=0
sprintf('2:条件平差BV+W=0')
n=input('请输入n=');
t=input('请输入t=');
B=input('请输入系数矩阵B=');
W=input('请输入自由项向量W=');
P=input('请输入权阵P=');
sprintf('开始计算')
K=-inv((B*inv(P)*B'))*W;
V=inv(P)*B'*K;
u=sqrt((V'*P*V)/(n-t));
Ql=inv(P)-inv(P)*B'*inv(B*inv(P)*B')*B*inv(P);
fid=fopen('data.txt','wt'); %写入文件路径，文件输出