基于主元素思想的Householder正交法解（矛盾）方程组

陆嵩 2015-12-06 02:10:53

1232

分类专栏： MATLAB与Mathematica 文章标签：计算机算法 householde

最后发布:2015-12-06 02:10:53首次发布:2015-12-06 02:10:53

本文链接：https://blog.csdn.net/lusongno1/article/details/50191591

版权

基于主元素思想的Householder正交法解（矛盾）方程组

[TOC]
[线性方程组]1求解已有比较成熟的算法，但仍然存在不少的问题亟待解决。比如一般的求解方法在求解百万阶的方程组时速度慢，由于计算机存储导致的舍入误差太大等问题，使人们想出了各种方法，另一方面也补充了LR分解和LDL $^{T}$ 以及QR分解在许多方面的局限性。

householder正交变换相比于高斯消元法和主元素方法，能够最大限度地保留原矩阵的性质，正交变换不改变矩阵的条件数
由于在正交变换化上三角的过程中，可能会由于对角线上产生零元素，导致Householder矩阵找不到，故借用列主元素的思想，我做了一个列方向上的排序
Householder方法进而可推广到矛盾方程组的求解，矛盾方程组无法求精确解，求得的是基于最小二乘上的一个近似解。

程序功能

设矩阵A为m行n列，求解复方程组AX=b，笔者已完成：

矩阵形状	满秩	不满秩	速度和精度
m=n	DONE	NOT YET	FAST
m>n	DONE	NOT YET	FAST
n>m	NOT YET	NOT YET	UNKNOWN

各程序块说明

主程序　
分各类情况求解各类方程组

基于主元素思想的householder正交法解（矛盾）方程组
clc;clear;
%% 情况一：nxn可逆复矩阵
n=5;
aug_A=randi(100,n,n+1)+randi(100,n,n+1)*i;%随机生成一个nxn的可逆复矩阵（增广)
x=OTH(aug_A)%调用household正交法函数解方程组
A=aug_A(1:n,1:n);
b=aug_A(1:n,n+1);
A\b%输出结果与内置方法作比较
%% 情况二：基于最小二乘思想解矛盾方程组（要求A满秩，不满秩的情况不讨论）
m=15;n=3;
aug_A=randi(100,m,n+1);
x=OTH(aug_A)
%% 情况三：m<n的多解方程组（包括A满秩和不满秩两种情况）
m=15;n=30;
aug_A=randi(100,m,n+1);
x=OTH(aug_A)
%% 情况四：nn矩阵A不可逆（无解或多解）
A=[3 1 2; 6 2 4; 3 3 5]
b=[1;6;3] 
aug_A=[A b]
OTH(aug_A)

OTH函数

Household正交化方法分解矩阵A

回代，逆向反解化简后的方程组

OTH.m
function x=OTH(aug_A)
[m,n_a]=size(aug_A);
n=n_a-1;
if m<n
    error('大兄弟，您输入的方程组个数不足啊！');
end
A=aug_A(1:m,1:n);
if m==n&&det(A)==0
    error('请输入可逆矩阵')
    return
end
b=aug_A(1:m,n_a);
for j=1:n-1
    M=[A(j:m,:) b(j:m)];
    M=sortrows(M,-j);
    A(j:m,:)=M(:,1:n);
    b(j:m)=M(:,n+1);
    p=householder(A(j:m,j));
    unimat=eye(j-1);
    p=blkdiag(unimat,p);
    A=p*A;
    b=p*b;
end
%% 回代求值
x=zeros(n,1);
for k=n:-1:1
    x(k)=(b(k)-sum(A(k,k+1:n)*x(k+1:n)))/A(k,k);
end
end

householder函数

输入一个任意的非零向量，返回一个Householder矩阵p，使得px=ke

1 $_{1}$ ,其中 e

1 $_{1}$ 表示输入向量的第一个分量的单位向量，一般取1。

Householder.m
function p=householder(x)
l=length(x);
x=reshape(x,l,1);
sigma=norm(x,2);
abs_x1=abs(x(1));
re_x1=real(x(1));
im_x1=imag(x(1));
alpha=atan(im_x1/re_x1);
if x(1)==0 
 alpha=0;
end
if re_x1<0
    if im_x1>=0
        alpha=alpha+pi;
    else
        alpha=alpha-pi;
    end
end
e1=[1;zeros(l-1,1)];
e_ialpha=exp(1i*alpha);(1+i)/sqrt(2);
k=-sigma*e_ialpha;
w=(x-k*e1)/sqrt(2*sigma^2+2*sigma*abs_x1);
p=eye(l)-2*(w*w');
end

用到的数学知识

定理： $\forall n\in\mathbb C \quad x\neq$ 0,一定可以找到Householder矩阵p=I-2ww $^{H}$ (w $\in\mathbb C$ , $\Vert$ w $\Vert$ $_2$ =1）,使得px=ke $_1$ 。
需要注意方程组求解的稳定性，对病态方程组，要尤为注意。
正交矩阵不改变矩阵的条件数1.

调用关系图:

导师宋士仓简介 宋士仓，男，1963年1月生。1982年7月郑州大学数学系数学专业毕业留校任教，期间1986.9-1989.5国防科技大学攻读硕士学位,1999.9-2002.7郑州大学读博士学位,2003.9-2005.8中国科学院数学与系统科学研究院做博士后研究。2004年晋升教授。主讲高等数学、计算方法、数学分析、复变函数、C语言、数据库等多门课程，是省级重点课程《数学分析》的主要成员。教学效果优秀，曾获郑州大学教学优秀青年奖（1996）、河南省职业技能大赛二等奖（2003）、郑州大学三育人积极分子等（2007，2012）、河南省优秀教师称号（2012）。参与原国家教委教改项目《高等数学CAI软件质疑与自测系统》的研制，作为主要参加者荣获河南省教学成果一等奖。目前从事的专业为计算数学专业，研究方向为微分方程数值解、有限元混合元方法、复合材料物理问题的多尺度方法.目前结合参与国家自然科学基金重点项目《空天飞行器材料与结构的性能评价及关键理论研究》和国家自然科学基金项目《各向异性有限元》的研究，发表了与此相关的十余篇文章