数值计算.docx
《数值计算.docx》由会员分享，可在线阅读，更多相关《数值计算.docx（42页珍藏版）》请在冰点文库上搜索。
数值计算
题目1
解方程
，其中
，要求
。
迭代原理：
此题涉及到复数方程组的求解，因此矩阵求逆时要采用复数矩阵的求逆方法,复矩阵求逆的全选主元高斯—约当法的算法同实矩阵的求逆算法，其中复矩阵分成实部与虚部两部分，即
算法中所有运算均为复数运算。
在求解时，首先要用牛顿迭代法（当然其他的迭代方法在可以求解的情况下也是可以的）求出一个近似解（不需要太精确），接下来用迭代改善法去求解，计算出残余向量
，再求解一次方程组
，改善：
令
，反复迭代改善，最终满足精度要求。
复数矩阵求逆的全选主元高斯约当法运算的过程:
可见，转化后的方程组是关于a,b,c,d四个未知量的方程组。
此时，转化后的方程组需化成严主对角占优的方程组，否则将不能使用迭代法进行运算.这里我们化简为:
流程图
运算结果及小结：
由上图中的计算结果可以看出x1=0.453448275862069+0.691379310344828i;
x2=0.406551724137931-0.171379310344828i.
（1）牛顿迭代法的改善法对于最终的结果是非常有效的，通过牛顿法求出一个不十分精确的值，经过改善后就可以得出最终的解。
（2）复数矩阵的运算先需要进行主对角占优，即按前面的方式进行计算即可。
通过这个方程的计算，熟悉了复数方程组的求解，理解了牛顿迭代的求解过程。
usingSystem;
usingSystem.Collections.Generic;
usingSystem.Text;
namespaceImprovementIteration
{
classProgram
{
staticvoidMain（string[]args）
{
//方程左边矩阵的实部矩阵
double[,]LR=newdouble[2,2];
//方程左边矩阵的虚部矩阵
double[,]LI=newdouble[2,2];
//方程右边矩阵的实部矩阵
double[]RR=newdouble[2];
//方程右边矩阵的虚部矩阵
double[]RI=newdouble[2];
//存储函数值
doubleFX1Imag=0;doubleFX2Imag=0;
doubleFX1Real=0;doubleFX2Real=0;
//存储初始变量的值
doubleX1R0=0;//实部doubleX2R0=0;
doubleX1I0=0;//虚部doubleX2I0=0;
//存储中间变量的值
doubleX1RK=0;doubleX2RK=0;
doubleX1IK=0;doubleX2IK=0;
doubleX1RK1=0;doubleX2RK1=0;
doubleX1IK1=0;doubleX2IK1=0;
//修正变量
doubleDX1RK=0;doubleDX2RK=0;
doubleDX1IK=0;doubleDX2IK=0;
//控制误差
doubleT=0;
for（inti=0;i<2;i++）
{
for（intj=0;j<2;j++）
{
LR[i,j]=0;LI[i,j]=0;
RR[i]=0;RI[i]=0;
}
}
//给矩阵赋值
LR[0,0]=8.0;LR[0,1]=1.0;LR[1,0]=-1.0;LR[1,1]=9.0;
LI[0,0]=1.0;LI[0,1]=-2;LI[1,0]=1.0;LI[1,1]=-3.0;
RR[0]=3.0;RR[1]=2.0;
RI[0]=5.0;RI[1]=-3.0;
//定义矩阵
MatrixLeftMatrixReal=newMatrix（2,2）;
MatrixLeftMatrixImag=newMatrix（2,2）;
MatrixRightMatrixReal=newMatrix（0,1）;
MatrixRightMatrixImag=newMatrix（0,1）;
//矩阵赋值
for（inti=0;i<2;i++）
{
for（intj=0;j<2;j++）
{
LeftMatrixImag.SetElement（i,j,LI[i,j]）;
LeftMatrixReal.SetElement（i,j,LR[i,j]）;
RightMatrixReal.SetElement（0,j,RR[j]）;
RightMatrixImag.SetElement（0,j,RI[j]）;
}
}
//计算复数方程组的雅克比矩阵的逆----程序见类Matrix附录B
LeftMatrixReal.InvertGaussJordan（LeftMatrixImag）;
//牛顿法求复数方程组的根
do
{//计算函数值f（x）
FX1Real=8.0*X1R0-1.0*X1I0+1.0*X2R0+2.0*X2I0-3;
FX1Imag=8.0*X1I0+1.0*X1R0+1.0*X2I0-2.0*X2R0-5;
FX2Real=-1.0*X1R0-1.0*X1I0+9.0*X2R0+3.0*X2I0-2;
FX2Imag=-1.0*X1I0+1.0*X1R0+9.0*X2I0-3.0*X2R0+3;
X1RK=X1R0-LeftMatrixReal.GetElement（0,0）*FX1Real-LeftMatrixImag.GetElement（0,0）*FX1Imag+LeftMatrixReal.GetElement（0,1）*FX2Real-LeftMatrixImag.GetElement（0,1）*FX2Imag;
X1IK=X1I0-LeftMatrixReal.GetElement（0,0）*FX1Imag+LeftMatrixImag.GetElement（0,0）*FX1Real+LeftMatrixReal.GetElement（0,1）*FX2Imag+LeftMatrixImag.GetElement（0,1）*FX2Real;
X2RK=X2R0-LeftMatrixReal.GetElement（1,0）*FX1Real-LeftMatrixImag.GetElement（1,0）*FX1Imag+LeftMatrixReal.GetElement（1,1）*FX2Real-LeftMatrixImag.GetElement（1,1）*FX2Imag;
X2IK=X2I0-LeftMatrixReal.GetElement（1,0）*FX1Imag+LeftMatrixImag.GetElement（1,0）*FX1Real+LeftMatrixReal.GetElement（1,1）*FX2Imag+LeftMatrixImag.GetElement（1,1）*FX2Real;
//计算残差r
RR[0]=3.0-（8.0*X1RK-1.0*X1IK+1.0*X2RK+2.0*X2IK）;
RR[1]=2.0-（-1.0*X1RK-1.0*X1IK+9.0*X2RK+3.0*X2IK）;
RI[0]=5.0-（8.0*X1IK+1.0*X1RK+1.0*X2IK-2.0*X2RK）;
RI[1]=-3.0-（-1.0*X1IK+1.0*X1RK+9.0*X2IK-3.0*X2RK）;
//计算修正量△x
DX1RK=LeftMatrixReal.GetElement（0,0）*RR[0]-LeftMatrixImag.GetElement（0,0）*RI[0]+LeftMatrixReal.GetElement（0,1）*RR[1]-LeftMatrixImag.GetElement（0,1）*RI[1];
DX1IK=LeftMatrixReal.GetElement（0,0）*RI[0]+LeftMatrixImag.GetElement（0,0）*RR[0]+LeftMatrixReal.GetElement（0,1）*RI[1]+LeftMatrixImag.GetElement（0,1）*RR[1];
DX2RK=LeftMatrixReal.GetElement（1,0）*RR[0]-LeftMatrixImag.GetElement（1,0）*RI[0]+LeftMatrixReal.GetElement（1,1）*RR[1]-LeftMatrixImag.GetElement（1,1）*RI[1];
DX2IK=LeftMatrixReal.GetElement（1,0）*RI[0]+LeftMatrixImag.GetElement（1,0）*RR[0]+LeftMatrixReal.GetElement（1,1）*RI[1]+LeftMatrixImag.GetElement（1,1）*RR[1];
//修正计算结果
X1RK1=X1RK+DX1RK;
X1IK1=X1IK+DX1IK;
X2RK1=X2RK+DX2RK;
X2IK1=X2IK+DX2IK;
T=X1RK1-X1R0+X1IK1-X1I0+X2RK1-X2R0+X2IK1-X2I0;
//结果输出
Console.WriteLine（"x1={0}+{1}i;x2={2}+{3}i",X1R0,X1I0,X2R0,X2I0）;
X1R0=X1RK1;
X1I0=X1IK1;
X2R0=X2RK1;
X2I0=X2IK1;
}while（Math.Abs（T）>0.000001）;
Console.Read（）;
}
}
}
运算结果：
x1=0+0i;x2=0+0i
x1=0.453448275862069+0.691379310344828i;x2=0.406551724137931+-0.171379310344828i
第二题
题目及要求：
其中
，要求
。
求x,y的值。
牛顿迭代法编程思路：
首先，先将方程改写为
的形式，
接下来，求
的雅克比矩阵，并将x的初值代入雅克比矩阵和
方程组构造的矩阵。
最后，解线性方程组
求出
，则
再计算出
如果
<10-5则满足精度要求，终止循环，如果不满足则重复第三步直至满足精度要求。
其中实矩阵的高斯约当法的算法原理为：
设A为n阶方阵，若存在n阶方阵B使AB=BA=I成立，则称A可逆，B是A的逆矩阵，记为A-1。
逆矩阵是惟一的。
高斯约当法（全选主元）求逆的步骤如下。
对于k从0到n-1做如下操作：
（1）从第k行、第k列开始的右下角子阵中选取绝对值最大的元素，并记为比元素所在的行号和列号，再通过行交换与列交换将它交换到主元素位置上（这一步成为全选主元）；
（2）
；
（3）
，j=0，1，2，...n-1;j
k;
（4）
（5）
根据在全选主元过程中所记录的行、列交换信息进行恢复，恢复的原则如下：
在全选主元的过程中，先交换的行、列后进行恢复；原来的行（列）交换用列（行）交换来恢复。
牛顿迭代法流程图
运算结果及小结：
由MicrosoftVisualstudio2008迭代出的结果可以看出当XK=1.62440573447633
YK=1.57622153939748时，可以满足精度要求。
小结：
此题花了我三天的时间，最后才解出来。
一开始我用一般迭代法，发现根号下的式子小于0。
在陈老师的指导下，我又变换了G（x）的形式，换了两种方法，结果都是G（x）>1，我以为是我计算出了问题，就编程序计算，发现都是G（X）>1，如下图所示：
最后我把初值换为X0=1.4,Y0=1.7时，便可以用一般迭代法求出结果，也能达到精度要求。
但是初值已经给定，我只能采用另外一种迭代方法，用牛顿迭代法不存在不收敛的问题，因此最终我选择用牛顿迭代法计算，并且最终在自己的不懈努力下得出了结果。
通过此题的编程，使我明白了一般迭代的局限性，不是所有的方程在所有的初值条件下都可以用一般迭代，在不能用一般迭代时可以考虑用牛顿迭代。
一般迭代的优点在于：
计算简便，不存在求逆运算，缺点在于：
同样一个方程，不同的初值可能是无解的。
牛顿迭代的优点在于：
所有的方程都可以求解。
缺点在于：
计算麻烦，存在求逆运算，占用较大存储空间。
附源程序：
usingSystem;
usingSystem.Collections.Generic;
usingSystem.Linq;
usingSystem.Text;
namespaceConsoleApplication20
{
classProgram
{
staticvoidMain（string[]args）
{//定义初值
doubleX0=1.2,Y0=2.4;
doubleXK=0,YK=0;
doubleXK1=0,YK1=0;
//存储函数值
double[]FXK=newdouble[2];
//修正变量
double[]D1=newdouble[2];
//定义迭代向量的一阶偏导数矩阵
double[,]A=newdouble[2,2];
doubleTAO=0;
//初始化迭代向量的一阶偏导数矩阵
for（inti=0;i<2;i++）
{
for（intj=0;j<2;j++）
{A[i,j]=0;}
}
for（inti=0;i<2;i++）
{FXK[i]=0;D1[i]=0;}
//求迭代向量的一阶偏导数矩阵
A[0,0]=2*X0+Y0*Y0*Y0;
A[0,1]=3*X0*Y0*Y0;
A[1,0]=6*X0;
A[1,1]=-3*Y0*Y0;
FXK[0]=X0*X0+X0*Y0*Y0*Y0-9;
FXK[1]=3*X0*X0-Y0*Y0*Y0-4;
XK=X0;
YK=Y0;
do
{A[0,0]=2*XK+YK*YK*YK;A[0,1]=3*XK*YK*YK;
A[1,0]=6*XK;A[1,1]=-3*YK*YK;
//求矩阵A的逆矩阵----程序见函数
InvertGaussJordan（A,2,2）;
FXK[0]=XK*XK+XK*YK*YK*YK-9;
FXK[1]=3*XK*XK-YK*YK*YK-4;
D1[0]=-（A[0,0]*FXK[0]+A[0,1]*FXK[1]）;
D1[1]=-（A[1,0]*FXK[0]+A[1,1]*FXK[1]）;
if（Math.Abs（D1[0]）>Math.Abs（D1[1]））
{TAO=Math.Abs（D1[0]）;}
else
{TAO=Math.Abs（D1[1]）;}
XK1=XK+D1[0];YK1=YK+D1[1];
XK=XK1;YK=YK1;
Console.WriteLine（"XK={0}",XK）;
Console.WriteLine（"YK={0}",YK）;
}while（Math.Abs（TAO）>1e-5）;
Console.Read（）;
}
计算结果：
XK=1.48861003861004
YK=1.7387727012727
XK=1.62063038602379
YK=1.58112711573076
XK=1.62440334180229
YK=1.57622794004076
XK=1.62440573447633
YK=1.57622153939748
题目三
，求
附近的特征值对应的特征向量，要求
迭代原理：
1、反幂法求解：
按照上述迭代公式，即可求出在特征值
附近的特征值以及对应的特征向量。
二、幂法求解：
k=1,2...
反幂法流程图
幂法流程图
运算结果及小结
（一）用反幂法得出的结果
（二）用幂法得到的结果
可以看出：
反幂法收敛速度较快，并且在同样的精度下，两种算法得出的特征值和特征向量也是不同的。
通过比较两种算法的结果和程序可以知道幂法和反幂法的特点：
幂法优点：
方法简单。
缺点在于：
它会把所有的特征值都求出来，增加了不必要的算法。
反幂法优点：
对于一个给定的近似特征值的方程，它是一个有效方法。
缺点在于：
求逆矩阵麻烦，占用存储空间较大。
源程序一（反幂法）
usingSystem;
usingSystem.Collections.Generic;
usingSystem.Linq;
usingSystem.Text;
namespaceConsoleApplication4
{
classProgram
{
staticvoidMain（string[]args）
{//定义矩阵
double[,]A=newdouble[3,3];
//定义初始向量
double[]V0=newdouble[3];
//定义迭代向量
double[]VK=newdouble[3];
double[]UK=newdouble[3];
//定义修正因子
doubleUMAX=0;doubleUMAX1=0;doubleT=0;//定义停机准则doubleCj=-13;
//初始化矩阵
for（inti=0;i<3;i++）
{
for（intj=0;j<3;j++）
{
A[i,j]=0;V0[i]=1;UK[i]=0;VK[i]=0;
}
}
//给矩阵赋值
A[0,0]=-12-Cj;A[0,1]=3;A[0,2]=3;
A[1,0]=3;A[1,1]=1-Cj;A[1,2]=-2;
A[2,0]=3;A[2,1]=-2;A[2,2]=7-Cj;
//求矩阵A的逆矩阵----程序见函数
InvertGaussJordan（A,3,3）;
do
{//计算Uk
for（inti=0;i<3;i++）
{
UK[i]=A[i,0]*V0[0]+A[i,1]*V0[1]+A[i,2]*V0[2];//
}
//求Uk的最大值
UMAX1=Math.Abs（UK[0]）;
for（inti=0;i<3;i++）
{
if（Math.Abs（UK[i]）>Math.Abs（UMAX1））
{UMAX1=UK[i];}
}
//求Vk
for（inti=0;i<3;i++）
{
VK[i]=UK[i]/UMAX1;
}
//计算两步迭代的差
T=UMAX1-UMAX;
//将K+1步的值赋给K--继续迭代
UMAX=UMAX1;
V0=VK;
//输出特征值和特征向量
Console.WriteLine（"λi={0}",Cj+1.0/UMAX）;
Console.WriteLine（"VK=（{0},{1},{2}）",VK[0],VK[1],VK[2]）;
}while（Math.Abs（T）>0.00001）;
Console.Read（）;
}
}
}
λi=-12.5925925925926
VK=（-1,0.271604938271605,0.197530864197531）
λi=-12.782470703125
VK=（1,-0.234537760416667,-0.171305338541667）
λi=-12.7797813584451
VK=（-1,0.235114344211104,0.171625202973873）
λi=-12.7798205952234
VK=（1,-0.235105350009682,-0.171621118249195）
λi=-12.7798200151753
VK=（-1,0.235105489426195,0.171621172182034）
λi=-12.7798200238357
VK=（1,-0.235105487274058,-0.171621171447389）
源程序二（幂法）
usingSystem;
usingSystem.Collections.Generic;
usingSystem.Linq;
usingSystem.Text;
namespaceConsoleApplication4
{
classProgram
{
staticvoidMain（string[]args）
{//定义矩阵
double[,]A=newdouble[3,3];
//定义初始向量
double[]V0=newdouble[3];
//定义迭代向量
double[]VK=newdouble[
升级会员 