《数值分析》上机实习报告.docx
- 文档编号:4311594
- 上传时间:2022-11-29
- 格式:DOCX
- 页数:19
- 大小:22.95KB
《数值分析》上机实习报告.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;j sr1+=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;k for(k=0;k for(l=0;l for(k=0;k for(k=0;k for(k=0;k for(l=0;l sr1=0;kr=0;ar=0.0; } system("cls"); for(k=0;k { for(l=0;l printf("\n"); } for(i=0;i for(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;m for(i=0;i for(j=0;j for(i=0;i for(j=0;j for(i=0;i for(j=0;j for(i=0;i for(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;i printf("%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;i for(j=0;j for(k=0;k { for(j=k+1;j for(i=k+1;i for(j=k+1;j for(i=k+1;i for(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;j h+=a2[i][j]*x[j]; x[i]=(b[i]-h)/a2[i][i]; h=0; } printf("\n高斯法结果: \n"); printf("\n"); for(i=0;i printf("%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;i for(i=0;i for(j=0;j for(i=0;i for(tt=0;tt { for(m=0;m for(k=0;k { for(j=k+1;j for(i=k+1;i for(j=k+1;j for(i=k+1;i for(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;j h+=a3[i][j]*yy[j]; yy[i]=(z[i]-h)/a3[i][i]; h=0; } for(m=0;m if(fabs(mk) for(m=0;m if(tt<(N-1))mk=0; } t=t1+1/mk; printf("\n反幂法的结果是: "); printf("\nd=%f\n",t);/*输出反幂法的解*/ for(m=0;m return0; } 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;i x[i+1]=1+i*h; b[0]=2*h*y[0]; b[N-1]=2*h*y[N-1]; for(i=1;i 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]; 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;i k[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;i sum=sum+c[i]*d[i]; printf("f(4.563)的值为: \n"); printf("f(4.536)=%10.8f\n",sum); for(i=0;i k2[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
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 数值分析 数值 分析 上机 实习 报告