《数值分析》上机实习报告.docx
《《数值分析》上机实习报告.docx》由会员分享，可在线阅读，更多相关《《数值分析》上机实习报告.docx（19页珍藏版）》请在冰豆网上搜索。
《数值分析》上机实习报告
数值分析程序设计作业
一.housholder变换.求方程的解
二.求特征值
一,程序要求:
（1）用Household变换，把A化为三对角阵B（并打印B）.
（2）用超松弛法求解Bx=b（取松弛因子ω=1.4,X（0）=0,迭代9次）
（3）用列主元消去法求解Bx=b。
（4）用对分法求B（即A）的小于23且最接近23的特征值的近似值（误差小于0.1）
（5）用反幂法求B的该特征值的更精确近似值及相应的特征向量.
二，解题算法
1：
Household法
（1）令A0=A,aij
（1）=aij，已知Ar-1即Ar-1=（aij（r））
（2）Sr=（（air（r））2）1/2
（3）αr=Sr2+|a（r）r+1,r|Sr
ur=[0,0,a（r）r+1,r+Sign（a（r）r+1,r）Sr,a（r）r+2,r,…,anr（r）]T
（4）yr=Ar-1ur/αr
（5）kr=urTyr/2αr
（6）qr=yr-krur
（7）Ar=Ar-1-（qrurT+urqrT）r=（1,2,……n-2）
2：
超松弛法
其基本思想是迭代，即在高斯方法已求出x（m）,x（m-1）的基础上，组合新的序列，从而加快收敛速度。
其算式是：
xi（m）=（1-ω）xi（m-1）+ω（bijxi（m）+xj（m-1）+gi）
其中ω是超松弛因子，当ω>1时，可以加快收敛速度。
3：
高斯消去法：
对矩阵作恰当的调整，选取绝对值尽量大的元素作为主元素。
然后把矩阵化
为上三角阵，再进行回代，求出方程的解。
对分法和反幂法求解特征值及特征向量
4.对分法:
对于对称的三对角方阵,其特征值必落在[-‖B‖,‖B‖]中,其中‖B‖是其任意范数.通过同号数的计算,可知在[0,‖B‖]中特征根的个数,在不断对分,就可得特征根所在的区间,如在[lk,uk]中有一根λk,且[lk,uk]很小,就可取lk+uk/2作为λk的近似值.
5.反幂法:
基于乘幂法中求A-1不太容易,因此不直接用A-1xk=xk+1,而使用Axk+1=xk.用求解线性方程组求出xk+1,再标准化,就可得到绝对值最小的特征值与相应的特征向量,这就是反代法（也称反幂法）.
三,程序清单
#include"math.h"
#include"io.h"
#include"stdio.h"
#include"string.h"
#defineN9
floatsign（floatx）;
voidguss（floata[N][N]）;
floatfanmi（floata[N][N]）;
floatsor（floata[N][N],floatb[N]）;
floatduifen（）;
voidhousholder（）;
intf（floatx）;
floate[N],o[N];
floatapprox;
floatb[N]={2.1874369,33.992318,-25.173417,0.84671695,1.784317,-86.612343,1.1101230,4.719345,-5.6784392};
floata[N][N]=
{
12.38412,2.115237,-1.061074,1.112336,-0.113584,0.718719,1.742382,3.067813,-2.031743,
2.115237,19.141823,-3.125432,-1.012345,2.189736,1.563849,-0.784165,1.112348,3.123124,
-1.061074,-3.125432,15.567914,3.123848,2.031454,1.836742,-1.056781,0.336993,-1.010103,
1.112336,-1.012345,3.123848,27.108437,4.101011,-3.741856,2.101023,-0.71828,-0.037585,
-0.113584,2.189736,2.031454,4.101011,19.897918,0.431637,-3.111223,2.121314,1.784317,
0.718719,1.563849,1.836742,-3.741856,0.431637,9.789365,-0.103458,-1.103456,0.238417,
1.742382,-0.784156,-1.056781,2.101023,-3.111223,-0.103458,14.713846,3.123789,-2.213474,
3.067813,1.112348,0.336993,-0.71828,2.121314,-1.103456,3.123789,30.719334,4.446782,
-2.031743,3.123124,-1.010103,-0.037585,1.784317,0.238417,-2.213474,4.446782,40.00001
};
voidhousholder（）
{
floatUr[N],Yr[N],qr[N];
inti,j,l,k;
floatsr1=0,sr=0.0,ar=0.0,*p,kr=0.0;
for（i=0;i{
for（j=i+1;jsr1+=a[j][i]*a[j][i];
sr=sqrt（sr1）;
ar=sr1+fabs（a[i+1][i]）*sr;
for（l=0;l{
if（l<=i）Ur[l]=0.0;
elseif（l==i+1）Ur[l]=a[l][i]+sign（a[l][i]）*sr;
elseUr[l]=a[l][i];
}
for（k=0;kfor（k=0;kfor（l=0;lfor（k=0;kfor（k=0;kfor（k=0;kfor（l=0;lsr1=0;kr=0;ar=0.0;
}
system（"cls"）;
for（k=0;k{
for（l=0;lprintf（"\n"）;
}
for（i=0;ifor（i=0;i}
floatsor（floata[N][N],floatb[N]）/*超松驰法求解的子函数*/
{
floata1[N][N];
inti,j,k,m;
floatx[N],temp[N][N];
floatw=1.4,h=0,g=0;
for（m=0;mfor（i=0;ifor（j=0;jfor（i=0;ifor（j=0;jfor（i=0;ifor（j=0;jfor（i=0;ifor（k=0;k{
for（i=1;i<=N;i++）
{
for（j=1;j<=N;j++）
{
if（j<=i-1）h+=temp[i-1][j-1]*x[j-1];
elseif（j==i）h+=0;
elseh+=temp[i-1][j-1]*x[j-1];
}
g=（h+b[i-1]）*w;
x[i-1]=（1-w）*x[i-1]+g;
h=0.0;g=0.0;
}
}
printf（"超松驰法结果：
:
\n"）;
printf（"\n"）;
for（i=0;iprintf（"%f,",x[i]）;/*输出超松驰法求解的特征向量*/
return0;
}
floatsign（floatx）/*子函数符号函数*/
{
if（x>0.0）return1.0;
elseif（x<0.0）return-1.0;
elsereturn0.0;
}
voidguss（floata[N][N]）/*高斯法求解的子函数*/
{
floata2[N][N];
inti,j,k;
floatl[N]={0,},x[N],h=0;
floatb[N]={2.1874369,33.992318,-25.173417,0.84671695,1.784317,-86.612343,1.1101230,4.719345,-5.6784392};
for（i=0;ifor（j=0;jfor（k=0;k{
for（j=k+1;jfor（i=k+1;ifor（j=k+1;jfor（i=k+1;ifor（i=k+1;i}
x[N-1]=b[N-1]/a2[N-1][N-1];
for（i=N-2;i>-1;i--）
{
for（j=i+1;jh+=a2[i][j]*x[j];
x[i]=（b[i]-h）/a2[i][i];
h=0;
}
printf（"\n高斯法结果:
\n"）;
printf（"\n"）;
for（i=0;iprintf（"%f,",x[i]）;/*输出高斯法求解的特征向量*/
printf（"\n"）;
return;
}
floatduifen（）/*对分法求取方程的解的子函数*/
{
inti,j;
floatss1=0,ss2=23,ss=0,h=0,g=0;
floatl,ll;
h=f（ss2）;
for（i=0;1;i++）
{
ss1=ss;
g=f（ss）;
ss=0.5*（ss+23）;
l=0.5*（ss1+23）;
if（g-h==1）break;
}
do{
if（f（ss）-h==0）{h=f（ss）;ss2=ss;}
elsess1=ss;
ss=0.5*（ss1+ss2）;
ll=0.5*（ss1+ss2）;
if（fabs（ll-l）<0.01）break;
elsel=ll;
}while
（1）;
returnll;
}
intf（floatx）/*求取每次新的对分点的函数值*/
{
inti,j,hh=0,n=9;
floatsa[9];
sa[0]=o[0]-x;
if（sa[0]==0）sa[1]=e[0]*e[0];
elsesa[1]=o[1]-x-e[0]*e[0]/sa[0];
for（i=2;i{
if（sa[i-1]*sa[i-2]!
=0）
sa[i]=o[i]-x-e[i-1]*e[i-1]/sa[i-1];
elseif（sa[i-2]==0）sa[i]=o[i]-x;
if（sa[i-1]==0）sa[i]=-1;
}
for（i=0;i{
if（sa[i]>=0）hh+=1;
elsehh+=0;
}
returnhh;
}
floatfanmi（floata[N][N]）/*反幂法求解方程的解的子函数*/
{
floatmk=0,t,t1;
floatz[N],yy[N],a3[N][N];
inti,j,k,m,tt;
floattemp[N][N];
floath=0;
floatl[N]={0,};
t1=duifen（）;/*调用对分法求取方程的解*/
for（i=0;ifor（i=0;ifor（j=0;jfor（i=0;ifor（tt=0;tt{
for（m=0;mfor（k=0;k{
for（j=k+1;jfor（i=k+1;ifor（j=k+1;jfor（i=k+1;ifor（i=k+1;i}
yy[N-1]=z[N-1]/a3[N-1][N-1];
for（i=N-1;i>-1;i--）
{
for（j=i+1;jh+=a3[i][j]*yy[j];
yy[i]=（z[i]-h）/a3[i][i];
h=0;
}
for（m=0;mif（fabs（mk）for（m=0;mif（tt<（N-1））mk=0;
}
t=t1+1/mk;
printf（"\n反幂法的结果是:
"）;
printf（"\nd=%f\n",t）;/*输出反幂法的解*/
for（m=0;mreturn0;
}
main（）
{
housholder（）;
sor（a,b）;/*调用超松驰法求解*/
guss（a）;/*调用高斯法求解*/
approx=duifen（a）;/*调用对分法求解*/
printf（"\n对分法的结果是:
"）;
printf（"\nd=%f",approx）;/*输出对分法近似解*/
fanmi（a）;/*调用反幂法求解*/
printf（"\n第一题完"）;
}_
四,运行结果
housholder变换为:
12.3841,-4.8931,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,
-4.8931,25.3984,6.4941,0.0000,-0.0000,0.0000,-0.0000,-0.0000,0.0000,
0.0000,6.4941,20.6115,8.2439,0.0000,0.0000,0.0000,0.0000,0.0000,
0.0000,0.0000,8.2439,23.4228,-13.8801,0.0000,0.0000,0.0000,0.0000,
0.0000,-0.0000,0.0000,-13.8801,29.6983,4.5345,0.0000,-0.0000,0.0000,
0.0000,0.0000,0.0000,0.0000,4.5345,16.0061,4.8814,0.0000,0.0000,
0.0000,0.0000,0.0000,0.0000,0.0000,4.8814,26.0134,-4.5036,0.0000,
0.0000,-0.0000,0.0000,0.0000,-0.0000,0.0000,-4.5036,21.2540,4.5045,
0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,4.5045,14.5341,
超松驰法的结果:
1.073410,2.272582,-2.856601,2.292513,2.112164,-6.422586,1.357799,0.634258,-0.587043,
高斯法的结果:
1.075800,2.275746,-2.855515,2.293098,2.112633,-6.423838,1.357920,0.634243,-0.587267,
对分法的结果:
d=21.916260
反幂法的结果:
d=21.928205
0.157064,-0.306357,0.282195,0.285922,0.198636,0.533750,0.462840,1.000000,0.611851,
五．问题讨论
问题讨论
1.SOR方法的矩阵形式为：
X（m）=（E-ωL）-1（（1-ω）E+ωR）x（m-1）+（E-ωL）-1ωg
若记Lω=（E-ωL）-1（（1-ω）E+ωR），SOR收敛的充要条件是S（Lω）<1.且若A为对称正定阵，则当松弛因子ω满足0<ω<2时，SOR方法收敛。
此题中矩阵B是对称正定阵，且是三对角的，所以选择合适的松驰因子ω，收敛速度是很快的。
2.对分法简单可靠,数值稳定性较高,对于求少量几个特征值特别适宜.
3.反幂法是结合对分法使用的,所以只要近似值较恰当,其收敛速度是惊人的.只要迭代两次就可得到较满意的结果.但在运用中需把一般矩阵化为Hessebberg阵,然后进行反迭代。
三.用三次样条插值求函数值
一,程序要求
已知十点函数值及起点和终点的导数值,要求用三次样条插值求f（4.563）及f'（3）的近似值.
二,程序算法
其基本思想是对均匀分划的插值函数的构造,三次样条函数空间中不取1,x,,x2,x3,（x-xj）+3为基函数,而取B样条函数Ω3（x-xj/h）为基函数.由于三次样条函数空间是N+3维的,故我们把分点扩大到X-1,XN+1,则任意三次样条函数可用Ω3（x-xj/h）线性组合来表示S（x）=cjΩ3（x-xj/h）这样对不同插值问题,若能确定cj由解的唯一性就能求得解S（x）.本题是属于问题一,由s（xi）=yi,I=1,2,…Ns'（x0）=y0',s'（xN）=yN'可得
S（xi）=cjΩ3（xi-xj/h）=yi
S'（x0）=1/hcjΩ3'（x0-xj/h）=y'0
S'（xN）=1/hcjΩ3'（xN-xj/h）=y'N
通过列主元消去法解出cj,在反代回去求出所需的函数值.
三,程序清单
#include
#include
#defineN12
floata[N][N]={{-2,-4,0,0,0,0,0,0,0,0,0,0},
{1,4,1,0,0,0,0,0,0,0,0,0},
{0,1,4,1,0,0,0,0,0,0,0,0},
{0,0,1,4,1,0,0,0,0,0,0,0},
{0,0,0,1,4,1,0,0,0,0,0,0},
{0,0,0,0,1,4,1,0,0,0,0,0},
{0,0,0,0,0,1,4,1,0,0,0,0},
{0,0,0,0,0,0,1,4,1,0,0,0},
{0,0,0,0,0,0,0,1,4,1,0,0},
{0,0,0,0,0,0,0,0,1,4,1,0},
{0,0,0,0,0,0,0,0,0,1,4,1},
{0,0,0,0,0,0,0,0,0,0,4,2}};
floatx[N],c[N],sum=0;
floaty[N]={1,0,0.69314718,1.0986123,1.3862944,1.6094378,1.7917595,1.9459101,2.079445,2.1972246,2.3025851,0.1};/*以上是预处理和定义全局变量*/
main（）
{inti;
floath=1.0,b[N],e[N],l[N],f[N],d[N],k[N],d1[N],d2[N],k1[N],k2[N],t,r;
x[1]=1;x[0]=x[1]-h;
for（i=1;ix[i+1]=1+i*h;
b[0]=2*h*y[0];
b[N-1]=2*h*y[N-1];
for（i=1;ib[i]=6*y[i];
b[0]=b[0]-b[1];b[N-1]=b[N-1]+b[N-2];
e[0]=a[0][0];f[0]=b[0];
for（i=0;i{l[i+1]=a[i+1][i]/e[i];
e[i+1]=a[i+1][i+1]-l[i+1]*a[i][i+1];
f[i+1]=b[i+1]-l[i+1]*f[i];}
c[N-1]=f[N-1]/e[N-1];
for（i=N-2;i>=0;i--）
c[i]=（f[i]-a[i][i+1]*c[i+1]）/e[i];
for（i=0;ik[i]=4.563-x[i];
for（i=0;i{if（fabs（k[i]）>=2）d[i]=0;
elseif（fabs（k[i]）<=1）
d[i]=（1/2.0）*pow（fabs（k[i]）,3）-k[i]*k[i]+（2/3.0）;
elseif（（fabs（k[i]）<2）&&（fabs（k[i]）>1））
d[i]=（-1/6.0）*pow（fabs（k[i]）,3）+k[i]*k[i]-2*fabs（k[i]）+（4/3.0）;}
for（i=0;isum=sum+c[i]*d[i];
printf（"f（4.563）的值为:
\n"）;
printf（"f（4.536）=%10.8f\n",sum）;
for（i=0;ik2[i]=3.5-x[i];
for（i=0;i{if（fabs（k2[i]）>=1.5）d2[i]=0;
elseif（fabs（k2[i]）<=（1/2.0））
d2[i]=-（k2[i]*k2[i]）+（3/4.0）;
elsei
升级会员 