《数值分析》上机实习报告Word格式.docx
- 文档编号:17233943
- 上传时间:2022-11-29
- 格式:DOCX
- 页数:19
- 大小:22.95KB
《数值分析》上机实习报告Word格式.docx
《《数值分析》上机实习报告Word格式.docx》由会员分享,可在线阅读,更多相关《《数值分析》上机实习报告Word格式.docx(19页珍藏版)》请在冰豆网上搜索。
5.反幂法:
基于乘幂法中求A-1不太容易,因此不直接用A-1xk=xk+1,而使用Axk+1=xk.用求解线性方程组求出xk+1,再标准化,就可得到绝对值最小的特征值与相应的特征向量,这就是反代法(也称反幂法).
三,程序清单
#include"
math.h"
io.h"
stdio.h"
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<
N-2;
i++)
for(j=i+1;
j<
N;
j++)
sr1+=a[j][i]*a[j][i];
sr=sqrt(sr1);
ar=sr1+fabs(a[i+1][i])*sr;
for(l=0;
l<
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;
k<
k++){Yr[k]=0;
k++)
for(l=0;
l++)Yr[k]+=a[k][l]/ar*Ur[l];
k++)kr+=Ur[k]*0.5/ar*Yr[k];
for(k=0;
k++)qr[k]=Yr[k]-kr*Ur[k];
l++)a[k][l]-=(qr[k]*Ur[l]+Ur[k]*qr[l]);
sr1=0;
kr=0;
ar=0.0;
system("
cls"
);
l++)printf("
%.4f,"
a[k][l]);
printf("
\n"
i++)o[i]=a[i][i];
N-1;
i++)e[i]=a[i][i+1];
}
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;
m<
m++)x[m]=0;
for(j=0;
j++)a1[i][j]=a[i][j];
j++)temp[i][j]=0.0;
j++)temp[i][j]=-a1[i][j]/a1[i][i];
i++)b[i]/=a1[i][i];
for(i=1;
=N;
for(j=1;
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;
超松驰法结果:
:
%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;
j++)a2[i][j]=a[i][j];
for(j=k+1;
j++)l[j]=a2[j][k]/a2[k][k];
for(i=k+1;
j++)a2[i][j]-=l[i]*a2[k][j];
i++)a2[i][k]=0;
i++)b[i]-=l[i]*b[k];
x[N-1]=b[N-1]/a2[N-1][N-1];
for(i=N-2;
i>
-1;
i--)
h+=a2[i][j]*x[j];
x[i]=(b[i]-h)/a2[i][i];
h=0;
\n高斯法结果:
/*输出高斯法求解的特征向量*/
return;
floatduifen()/*对分法求取方程的解的子函数*/
inti,j;
floatss1=0,ss2=23,ss=0,h=0,g=0;
floatl,ll;
h=f(ss2);
1;
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;
n;
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;
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();
/*调用对分法求取方程的解*/
i++)z[i]=1;
j++)a3[i][j]=a[i][j];
i++)a3[i][i]-=t1;
for(tt=0;
tt<
tt++)
m++)yy[m]=0;
j++)l[j]=a3[j][k]/a3[k][k];
j++)a3[i][j]-=l[i]*a3[k][j];
i++)a3[i][k]=0;
i++)z[i]-=l[i]*z[k];
yy[N-1]=z[N-1]/a3[N-1][N-1];
for(i=N-1;
h+=a3[i][j]*yy[j];
yy[i]=(z[i]-h)/a3[i][i];
m++)
if(fabs(mk)<
fabs(yy[m]))mk=yy[m];
m++)z[m]=yy[m]/mk;
if(tt<
(N-1))mk=0;
t=t1+1/mk;
\n反幂法的结果是:
"
\nd=%f\n"
t);
/*输出反幂法的解*/
m++)printf("
z[m]);
/*输出反幂法的解的特征向量*/
main()
{
housholder();
sor(a,b);
/*调用超松驰法求解*/
guss(a);
/*调用高斯法求解*/
approx=duifen(a);
/*调用对分法求解*/
\n对分法的结果是:
\nd=%f"
approx);
/*输出对分法近似解*/
fanmi(a);
/*调用反幂法求解*/
\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'
(xN)=1/hcjΩ3'
(xN-xj/h)=y'
N
通过列主元消去法解出cj,在反代回去求出所需的函数值.
#include<
stdio.h>
math.h>
#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;
x[i+1]=1+i*h;
b[0]=2*h*y[0];
b[N-1]=2*h*y[N-1];
b[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];
{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];
=0;
c[i]=(f[i]-a[i][i+1]*c[i+1])/e[i];
k[i]=4.563-x[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);
sum=sum+c[i]*d[i];
f(4.563)的值为:
f(4.536)=%10.8f\n"
sum);
k2[i]=3.5-x[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
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 数值分析 数值 分析 上机 实习 报告