北航研究生数值分析编程大作业2.docx
- 文档编号:2631615
- 上传时间:2023-05-04
- 格式:DOCX
- 页数:40
- 大小:1.74MB
北航研究生数值分析编程大作业2.docx
《北航研究生数值分析编程大作业2.docx》由会员分享,可在线阅读,更多相关《北航研究生数值分析编程大作业2.docx(40页珍藏版)》请在冰点文库上搜索。
北航研究生数值分析编程大作业2
数值分析
计算实习二
姓名:
徐航
学号:
SY1105517
学院:
航空科学与工程学院
具体源程序见光盘
一、开发平台的选择
本题中选择MicrosoftVisualStudio2008作为开发平台,利用VC++作为开发环境,编了一个基于MFC的程序。
可以使用VS2008打开“源文件”文件夹中的“QR_2.sln”查看或运行源文件。
本程序不仅仅是界面不同,编程的思想和运算的速度都有提高。
其优点是通过设计类,将所求数组、维数等封装为成员变量。
将矩阵的加法、乘法、数乘等运算符进行重载,使得在遇到两个矩阵相乘时只需写乘号而无需编写循环代码。
将带双步位移的QR分解、矩阵的拟上三角化等功能封装成类的成员函数,在定义了对象后可以任意调用。
这都是利用了封装性。
在定义类的对象之后可以方便地调用这些成员变量和成员函数。
本程序可以很容易地升级为通用性的程序,用户利用可执行文件只需要输入几个参数就可以运算同一类的任意方程,而用户不用像C语言那样改动源代码。
应用面向对象的编程思想还有一个最大的优点是可扩展性。
想添加或改变功能只需要改动成员函数就行了而不会影响到程序的其他部分。
同时,利用MFC还可以开发出很友好的程序界面。
本次程序的初始化界面如下:
图1程序初始化界面
初始化的界面中会显示原矩阵A。
二、算法的设计方案
1、矩阵的存储方案
本方案采用10*10的矩阵来存储对象,虽然在运算时下标总要减1(由于C++中数组的下标是从0开始的),但出于节省内存空间的考虑,还是没有采用11*11的简便方法。
2、求解的总体思路
1)界面进行初始化之后点击
按钮。
激发Windows的消息响应。
在消息响应函数中先判断是否已经计算过了。
若没有计算过则进行计算。
否则弹出对话框如图:
图2当重复计算时弹出的对话框
此时若单击按钮
则会因为还没有进行计算而弹出对话框:
图3当还没有进行计算就想显示特征值时弹出的对话框
当单击界面上的按钮
则会执行退出操作。
2)计算过程中首先对矩阵A进行拟上三角化,本程序中是调用子函数Get_Hessenberg()。
之后得到拟上三角阵A(n-1),并将它显示在控件中。
3)之后对矩阵A进行QR分解。
这里QR分解用到的算法是教材3.3.1节中的算法。
这里调用函数Pure_QR()完成了这一功能的。
函数执行完将得到Q矩阵和R矩阵。
利用重载的运算符“*”进行矩阵与矩阵的乘法“R*Q”将得到RQ矩阵。
将Q、R和RQ这三个矩阵显示在控件中。
4)调用DoubleQR()成员函数对拟上三角矩阵A(n-1)进行带双步位移的QR分解。
程序在无穷循环中应用switch语句达到在步骤之间转换的目的。
其基本步骤如下:
(1)定义
,给定精度水平
和迭代最大次数L。
(2)令
。
K表示迭代的次数,m表示当前迭代时矩阵的阶数。
(3)如果
,则得到A的一个特征值
,置
(降阶),转(4);否则转(5)。
(4)如果
或2则得到A的最后一个或两个特征值,转(9);否则转(3)。
(5)如果
则得到A的两个特征值
和
,置
(降阶),转(4);否则转(6)。
(6)如果k=L,则计算终止,未得到A的全部特征值;否则转(7)。
这一步以为着迭代达到了上限,算法失效。
(7)记
计算
(8)置
,转(3)。
(9)A的全部特征值已计算完毕,停止计算。
将得到的迭代结构矩阵——分块上三角阵显示在控件中。
5)定义子窗体的一个对象,调用MFC中的DoModal()函数显示出子窗体并将所有特征值显示在子窗体初始化的控件中。
6)调用Get_Eigenvector()函数,利用列主元素的高斯消去法求解特征方程,得到每一个实特征值所对应的特征向量。
并将这些特征向量显示在控件中。
此时显示特征值、特征向量的窗体为:
图4显示特征值和特征向量的窗体
三、源代码
在光盘中有全部的源代码,请用MicrosoftVisualStudio2008打开。
同时还有一个课执行文件可以直接执行。
由于有些程序段是MFC直接生成的,这里没有把那部分程序打出。
1、Array类的设计
类的头文件为:
#pragmaonce//使头文件只被编译一次
/*常量ROW_N表示矩阵A行的维数*/
#defineROW_N10
/*常量COL_N表示矩阵A列的维数*/
#defineCOL_N10
//定义精度
#defineACCURACY1.0e-12
classArray//存放所求矩阵A的类。
类的成员函数中定义了对类的基本操作。
{
public:
//声明构造函数和析构函数
Array(void);//默认的构造函数
Array(intm,intn);//重载构造函数
~Array(void);//默认的析构函数
public:
//声明成员变量
//矩阵的行数
intm_Row;
//矩阵的列数
intm_Col;
//A存储矩阵
doublem_Array[ROW_N][COL_N];
//测试
doubletest;
public:
//声明成员函数
//重载加法运算符,定义矩阵与矩阵的加法运算
Arrayoperator+(Array&array_b);
//重载减法运算符,定义矩阵与矩阵的乘法运算
Arrayoperator-(Array&array_b);
//重载乘法运算符,定义矩阵与矩阵的乘法运算
Arrayoperator*(Array&array_b);
//重载乘法运算符,定义矩阵的数乘运算
Arrayoperator*(doublenum);
//符号函数
doublesgn(doublenum);
//判断某一列从第r个元素开始以下的元素是否为零
boolZero_or_Not(inti,intr);
//归零函数,将小于精度的数设为零
voidSet_Zero();
//计算矩阵的转置
ArrayArray_T();
//计算向量中一部分元素组成的新向量的二范数
doubleGet_Norm(intr,inti1,inti2);
//将一维数组转化为数
doubleChange_Array_to_Num();
//将原矩阵A进行拟上三角化
voidGet_Hessenberg();
//将向量重新初始化
voidRecover_Vector(Array&array_b);
//进行QR分解的函数
doubleQR(intm,Array&C);
//双步位移的QR分解
intDoubleQR(Array&lambda);
//进行纯QR分解的函数
doublePure_QR(Array&C);
//找出按列主元素
intGet_max_in_line(Array&A,intk);
//交换矩阵的两行
intExchange_Line(Array&A,inta,intb,intk);
//求矩阵的特征向量
intGet_Eigenvector(doublelambda,Array&ev);
};
类的实现文件为:
#include"StdAfx.h"
#include"Array.h"
#include
#include"QR_2Dlg.h"
Array:
:
Array(void)//默认的构造函数
{
}
Array:
:
Array(intm,intn)
{
m_Row=m;
m_Col=n;
test=100.0;//测试用的变量
}
Array:
:
~Array(void)
{
}
//重载加法运算符,定义矩阵与矩阵的加法运算
ArrayArray:
:
operator+(Array&array_b)//array_b为“+”后面的数组
{
inti,j;
Arraytemp(this->m_Row,this->m_Col);//用来存储结果的数组
for(i=0;i for(j=0;j temp.m_Array[i][j]=this->m_Array[i][j]+array_b.m_Array[i][j]; returntemp; } //重载乘法运算符,定义矩阵与矩阵的乘法运算 ArrayArray: : operator*(Array&array_b)//array_b为“*”后面的对象 { inti,j,n; Arraytemp(this->m_Row,array_b.m_Col);//temp为结果数组 for(i=0;i for(j=0;j { temp.m_Array[i][j]=0.0; for(n=0;n temp.m_Array[i][j]+=this->m_Array[i][n]*array_b.m_Array[n][j]; } returntemp; } //重载乘法运算符,定义矩阵的数乘运算 ArrayArray: : operator*(doublenum) { inti,j; Arraytemp(this->m_Row,this->m_Col); for(i=0;i for(j=0;j temp.m_Array[i][j]=num*m_Array[i][j]; returntemp; } //重载减法运算符,定义矩阵与矩阵的乘法运算 ArrayArray: : operator-(Array&array_b) { return(*this)+(array_b*(-1)); } //符号函数 doubleArray: : sgn(doublenum) { if(fabs(num)>ACCURACY) returnnum>0.0? 1.0: (-1.0); else return0.0; } //判断某一列从第r个元素开始以下的元素是否为零,只要有一个不为零就返回假 boolArray: : Zero_or_Not(inti,intr) { intj; for(j=i;j<=m_Row;j++) if(fabs(m_Array[j-1][r-1])>1e-12)returnfalse; returntrue; } //归零函数,将小于精度的数设为零 voidArray: : Set_Zero() { for(inti=0;i for(intj=0;j if(fabs(m_Array[i][j]) this->m_Array[i][j]=0.0; } //计算矩阵的转置 ArrayArray: : Array_T() { Arraytemp(this->m_Col,this->m_Row); for(inti=0;i for(intj=0;j temp.m_Array[i][j]=this->m_Array[j][i]; returntemp; } //计算向量中一部分元素组成的新向量的二范数 doubleArray: : Get_Norm(intr,inti1,inti2)//正确 { doubleanswer=0.0; for(inti=i1;i<=i2;i++) answer=answer+m_Array[i-1][r-1]*m_Array[i-1][r-1]; returnsqrt(answer); } //此方法用来获得一阶方阵的唯一元素的值 doubleArray: : Change_Array_to_Num() { returnm_Array[0][0]; } voidArray: : Recover_Vector(Array&array_b)//将向量中的每个元素赋值为0 { for(inti=0;i for(intj=0;j array_b.m_Array[i][j]=0.0; } //将原矩阵A进行拟上三角化 voidArray: : Get_Hessenberg() { intr=0,j=0; Arrayur(m_Col,1),pr(m_Col,1),qr(m_Col,1),wr(m_Col,1),tr(1,1); doubledr=0.0,cr=0.0,hr=0.0; for(r=1;r<=COL_N-2;r++) { if(Zero_or_Not(r+2,r)) continue; else { Recover_Vector(ur); dr=Get_Norm(r,r+1,m_Row); cr=-sgn(m_Array[r][r-1])*dr; hr=cr*(cr-m_Array[r][r-1]); ur.m_Array[r][0]=m_Array[r][r-1]-cr; for(inti=r+2;i<=m_Row;i++) ur.m_Array[i-1][0]=m_Array[i-1][r-1]; pr=(this->Array_T())*ur*(1/hr); qr=(*this)*ur*(1/hr); tr=pr.Array_T()*ur*(1/hr); wr=qr-ur*(tr.Change_Array_to_Num()); (*this)=(*this)-(wr*ur.Array_T())-(ur*pr.Array_T()); } } Set_Zero(); } //按照教材3.3.1节的步骤进行纯QR分解的函数 doubleArray: : Pure_QR(Array&C)//此时C存储的是Q而不是Mk { intr,j; doubledr=0.0,cr,hr; Arrayur(m_Row,1),pr(m_Row,1),wr(m_Row,1); for(r=1;r<=m_Row-1;r++) { if(Zero_or_Not(r+1,r)) continue; else { dr=0; for(j=r;j<=m_Row;j++) dr=sqrt(dr*dr+m_Array[j-1][r-1]*m_Array[j-1][r-1]); cr=-sgn(m_Array[r-1][r-1])*dr; hr=cr*(cr-m_Array[r-1][r-1]); for(inti=1;i<=m_Row;i++) ur.m_Array[i-1][0]=(i 0: m_Array[i-1][r-1]; ur.m_Array[r-1][0]=m_Array[r-1][r-1]-cr; wr=C*ur; C=C-wr*ur.Array_T()*(1/hr); pr=this->Array_T()*ur*(1/hr); (*this)=(*this)-ur*pr.Array_T(); } } return0; } //进行QR分解的函数(这是带双步位移的QR分解的子程序,和书上.3.1节中的纯QR分解的步骤有区别) doubleArray: : QR(intm,Array&C)//C存储的是Mk,原矩阵存储的是R { intr,j; doubledr=0.0,cr,hr; Arrayur(m,1),vr(m,1),pr(m,1),qr(m,1),wr(m,1),tr(1,1); for(r=1;r<=m-1;r++) { if(Zero_or_Not(r+1,r)) continue; else { dr=0; for(j=r;j<=m;j++) dr=sqrt(dr*dr+m_Array[j-1][r-1]*m_Array[j-1][r-1]); cr=-sgn(m_Array[r-1][r-1])*dr; hr=cr*(cr-m_Array[r-1][r-1]); for(inti=1;i<=m;i++) ur.m_Array[i-1][0]=(i 0: m_Array[i-1][r-1]; ur.m_Array[r-1][0]=m_Array[r-1][r-1]-cr; vr=this->Array_T()*ur*(1/hr); (*this)=(*this)-ur*vr.Array_T(); pr=C.Array_T()*ur*(1/hr); qr=C*ur*(1/hr); tr=pr.Array_T()*ur*(1/hr); wr=qr-ur*(tr.Change_Array_to_Num()); C=C-wr*ur.Array_T()-ur*pr.Array_T(); } } return0; } intArray: : DoubleQR(Array&eigenvalue)//带双步位移的QR分解 { intL=1000;//定义最大迭代次数 intk=1;//定义当前迭代次数 intm=m_Col;//第一次迭代时Ak的阶数是m_Col intn=3;//指示应该转到第几步 inti=1;//表示求到第几个特征值了 doubledelta=1;//解二元一次方程组时用到的判别数 while (1) { switch(n)//通过无穷循环中switch语句达到在步骤之间转换的目的 { case3: //第三步,若最后一行的倒数第二个数为零,则得到一个特征值 { if(fabs(m_Array[m-1][m-2])<=ACCURACY) { eigenvalue.m_Array[i-1][0]=m_Array[m-1][m-1];//存储特征值的实部 eigenvalue.m_Array[i-1][1]=0;//存储特征值的虚部 i++; m=m-1;//矩阵降一阶 n=4;//转第四步 break; } else { n=5;//否则转第五步 break; } } case4: //如果m=1或m=2,则得到最后两个特征值,转(9),否则转(3) { if(m==1) { eigenvalue.m_Array[i-1][0]=m_Array[m-1][m-1];//存储特征值的实部 eigenvalue.m_Array[i-1][1]=0;//存储特征值的虚部 n=9;//转() break; } elseif(m==2) { delta=(m_Array[m-2][m-2]+m_Array[m-1][m-1])*(m_Array[m-2][m-2]+m_Array[m-1][m-1]) -4*(m_Array[m-2][m-2]*m_Array[m-1][m-1]-m_Array[m-1][m-2]*m_Array[m-2][m-1]);//delta=b^2-4ac(在下面的运算中再开方) if(delta<0) { eigenvalue.m_Array[i-1][0]=(m_Array[m-2][m-2]+m_Array[m-1][m-1])/2.0; eigenvalue.m_Array[i-1][1]=sqrt(-delta)/2.0; i++; eigenvalue.m_Array[i-1][0]=(m_Array[m-2][m-2]+m_Array[m-1][m-1])/2.0; eigenvalue.m_Array[i-1][1]=-(sqrt(-delta))/2.0;i++; } else { eigenvalue.m_Array[i-1][0]=(m_Array[m-2][m-2]+m_Array[m-1][m-1]+sqrt(delta))/2.0; eigenvalue.m_Array[i-1][1]=0;i++; eigenvalue.m_Array[i-1][0]=(m_Array[m-2][m-2]+m_Array[m-1][m-1]-sqrt(delta))/2.0; eigenvalue.m_Array[i-1][1]=0;i++; } n=9;//转(9) break; } else { n=3; break; } } case5: //如果m-1行的第m-2个数为零,则得到两个特征值,转(4)。 否则转(6)。 { if(fabs(m_Array[m-2][m-3])<=ACCURACY) { delta=(m_Array[m-2][m-2]+m_Array[m-1][m-1])*(m_Array[m-2][m-2]+m_Array[m-1][m-1]) -4*(m_Array[m-2][m-2]*m_Array[m-1][m-1]-m_Array[m-1][m-2]*m_Array[m-2][m-1]); if(delta<0) { eigenvalue.m_Array[i-1][0]=(m_Array[m-2][m-2]+m_Array[m-1][m-1])/2.0;eigenvalue.m_Array[i-1][1]=sqrt(-delta)/2.0;i++; eigenvalue.m_Array[i-1][0]=(m_Array[m-2][m-2]+m_Array[m-1][m-1])/2.0;eigenvalue.m_Array[i-1][1]=
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 北航 研究生 数值 分析 编程 作业