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

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

目录

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

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

程序功能

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

矩阵形状满秩不满秩速度和精度
m=nDONENOT YETFAST
m>nDONENOT YETFAST
n>mNOT YETNOT YETUNKNOWN

各程序块说明

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

基于主元素思想的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 ,其中 e 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

用到的数学知识

  • 定理: nCx 0,一定可以找到Householder矩阵p=I-2ww H (w C , w 2 =1),使得px=ke 1
  • 需要注意方程组求解的稳定性,对病态方程组,要尤为注意。

  • 正交矩阵不改变矩阵的条件数1.

调用关系图:

Created with Raphaël 2.1.0 主函数 主函数 OTH.m OTH.m Householder.m Householder.m 嘿,小O儿,我要调用你。 小O愣了一下。 哥,老板说要开工了,休息一会开工。

导师宋士仓简介 宋士仓,男,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. 之前已说过,对于一些专业性术语,不做解释,请自行查找。
已标记关键词 清除标记
©️2020 CSDN 皮肤主题: 点我我会动 设计师:上身试试 返回首页