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

目录

[TOC]
[线性方程组]1求解已有比较成熟的算法，但仍然存在不少的问题亟待解决。比如一般的求解方法在求解百万阶的方程组时速度慢，由于计算机存储导致的舍入误差太大等问题，使人们想出了各种方法，另一方面也补充了LR分解和LDLT以及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.mfunction x=OTH(aug_A)[m,n_a]=size(aug_A);n=n_a-1;if m<n    error('大兄弟，您输入的方程组个数不足啊！');endA=aug_A(1:m,1:n);if m==n&&det(A)==0    error('请输入可逆矩阵')    returnendb=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);endend

householder函数

输入一个任意的非零向量，返回一个Householder矩阵p，使得px=ke1,其中 e1表示输入向量的第一个分量的单位向量，一般取1。

Householder.mfunction 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;endif re_x1<0    if im_x1>=0        alpha=alpha+pi;    else        alpha=alpha-pi;    endende1=[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

用到的数学知识

• 定理： ∀n∈Cx≠0,一定可以找到Householder矩阵p=I-2wwH(w∈C,∥w∥2=1）,使得px=ke1。

• 需要注意方程组求解的稳定性，对病态方程组，要尤为注意。

• 正交矩阵不改变矩阵的条件数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软件质疑与自测系统》的研制，作为主要参加者荣获河南省教学成果一等奖。 目前从事的专业为计算数学专业，研究方向为微分方程数值解、有限元混合元方法、复合材料物理问题的多尺度方法.目前结合参与国家自然科学基金重点项目《空天飞行器材料与结构的性能评价及关键理论研究》和国家自然科学基金项目《各向异性有限元》的研究，发表了与此相关的十余篇文章

---

---

1. 之前已说过，对于一些专业性术语，不做解释，请自行查找。 ↩