数值分析计算实习作业一.docx_第1页
数值分析计算实习作业一.docx_第2页
数值分析计算实习作业一.docx_第3页
数值分析计算实习作业一.docx_第4页
数值分析计算实习作业一.docx_第5页
已阅读5页,还剩6页未读 继续免费阅读

下载本文档

版权说明:本文档由用户提供并上传,收益归属内容提供方,若内容存在侵权,请进行举报或认领

文档简介

数 值 分 析计算实习题一学号:姓名: 院系: 2015年11月5日数值分析计算实习作业一一、分析1.1算法分析题目要求求出:1)特征值从小到大排列的最小特征值和最大特征值。2)特征值中模最小的特征值。3)靠近一组数的一组特征值。4)矩阵A的条件数cond(A)2。5)行列式detA。解决方法:1)若将所有行列式按模的大小排列则模最大的特征值一定是和中的一个,因此利用幂法求出模最大的特征值。然后利用带原点平移的幂法,将系数矩阵变为即将所有特征值都减去,则特征值按大小顺序排列的次序不变,模最大的特征值依然在整个排列的两端,再用一次幂法得到模最大的特征值,其中为带原点平移的幂法求出的特征值,最后两个特征值、比较大小,大的为,小的为。2)因为为按模最小的特征值,因此用反幂法可求的其特征值。3)因为靠近数,因此一定是所有的中模最小的,因此可利用带原点平移的反幂法求出特征值,此时的系数矩阵变为。4)条件数cond(A)2为模最小的特征值与模最大的特征值的比的绝对值,因此利用1和2中求出的和可解出条件数。5)可对矩阵A进行LU分解,即则,又因为矩阵L对角线元素为1,则=1,所以,U为上三角阵,行列式为对角线元素的乘积,因此可得A的行列式。1.2程序分析1.2.1 因为A为拟三角阵,储存时零元素不储存,因此将矩阵A压缩为5*501的矩阵CA的带内元素=C中的元素程序中A5501即为压缩后的矩阵。1.1.2 程序中的B5501为过渡矩阵,在幂法迭代、反幂法迭代以及LU分解中均用矩阵B来计算,计算之间对B进行适当的赋值。先令B=A,调用幂法函数求;再令B=A-I,调用幂法函数求出;再令B=A,调用反幂法函数求;在令B=,调用反幂法函数求,最后令B=A,将B进行LU分解,计算A的行列式。1.2.3幂法求解过程:1)取初始向量u0=1,1,1;2)进行k次迭代,k=1,2,33)对每一次迭代,计算上一次迭代中的的模(2范数),将单位化后赋值给,即,计算矩阵B与向量的乘积,计算与的内积。4)如果两次相邻两次迭代的满足:(允许误差),则结束迭代,的值可认为是矩阵B对应的模最大的特征值。如果不满足误差条件,则重复3)的迭代直到达到误差允许值。1.2.4 反幂法求解过程:1)取初始向量u0=1,1,1;对矩阵B进行LU分解。2)进行k次迭代,k=1,2,33)对每一次迭代,计算上一次迭代中的的模(2范数),将单位化后赋值给,即,利用LU分解后的矩阵B,求线性方程组,计算与的内积。4)如果两次相邻两次迭代的满足:(允许误差),则结束迭代,的值可认为是矩阵B对应的模最大的特征值。如果不满足误差条件,则重复3)的迭代直到达到误差允许值。1.2.5 拟三角矩阵的LU分解过程可参照数值分析书26页的算法得到。二、程序整个数值分析在c语言中的程序如下:#include #include double mifa();double fanmifa();void fenjie();/声明三个函数,mifa是用幂法求模最大的特征值,fanmifa是用反幂法求模最小的特征值,fenjie是对系数矩阵做LU分解double B5501=0;double A5501=0;/定义两个系数矩阵为全局变量,B作为求解过程中的系数矩阵使用double sigma=1e-12;/定义误差允许值/*主程序*/void main() int i,j,t; double lamk,lamm,min,max,miu,lam,x,lams; double cond2; double det=1; for(i=0;i=500;i+)/输入矩阵A j=i+1;x=0.1/j;A2i=(1.64-0.024*j)*sin(0.2*j)-0.64*exp(x);A3i=0.16;A4i=-0.064; for(i=1;i=500;i+) A1i=0.16; for(i=2;i=500;i+) A0i=-0.064; A3500=0; A4500=0; A4499=0; for(i=0;i=4;i+)/先使B=A for(j=0;j=500;j+) Bij=Aij; lamk=mifa();/幂法求模最大的特征值 for(i=0;i=4;i+)/赋值计算B=(A-lamk*I) for(j=0;j=500;j+) if(i=2) Bij=Aij-lamk; else Bij=Aij; lamm=mifa()+lamk;/带原点位移的幂法求位移后模最大的特征值 if(lamklamm) min=lamk;max=lamm; else min=lamm;max=lamk; printf(lam1=%13.11len,min);/显示12位有效数字 printf(lam501=%13.11len,max); for(i=0;i=4;i+)/赋值B=A for(j=0;j=500;j+) Bij=Aij; lams=fanmifa(); printf(lams=%13.11len,lams);/反幂法求模最小的特征值 for(t=1;t=39;t+)/循环求与miu接近的特征值 miu=min+t*(max-min)/40;for(i=0;i=4;i+)/赋值B=(A-miu*I) for(j=0;j=500;j+) if(i=2) Bij=Aij-miu; else Bij=Aij; lam=fanmifa()+miu;/用带原点的反幂法求与miu接近的特征值printf(lami%d=%13.11len,t,lam); cond2=fabs(lams/lamk);/计算A的条件数 printf(condA=%13.11len,cond2); for(i=0;i=4;i+)/赋值B=A for(j=0;j=500;j+) Bij=Aij; fenjie();/对B进行LU分解,其第三行的元素即为矩阵U的对角线元素,即乘积为A的行列式的值 for(i=0;i=500;i+) det=det*B2i; printf(det(A)=%13.11len,det);/计算输出A的行列式/*幂法程序*/double mifa()/幂法的函数int i,j,p,s; double u501,y501=0; double ita,a,beta1=0,beta2=0,b,c,lam; for(i=0;i=500;i+) ui=1; while(1)/开始迭代 a=0;for(i=0;i=500;i+) a=a+ui*ui;ita=sqrt(a);/计算2范数for(i=0;i=500;i+)/计算向量y(k-1) yi=ui/ita;for(i=0;i=500;i+) ui=0;for(i=2;i=502;i+)/计算u(k) if(i4) p=i; else p=4; if(i500) s=i; else s=500; for(j=i-p;j=s;j+) ui-2=ui-2+Bi-jj*yj; for(i=0;i=500;i+) beta2=beta2+yi*ui;/计算beta(T)b=fabs(beta2-beta1);c=fabs(beta2);if(b/c)=sigma)/如果达到精度则结束并返回值 lam=beta2; goto finish;else beta1=beta2; beta2=0; finish:return(lam);/*反幂法程序*/double fanmifa()/反幂法函数 int i,k,t,p; double u501,y501=0,d501=0; double ita,a,beta1=0,beta2=0,b,c,lam; for(i=0;i=500;i+) ui=1; fenjie();/对系数矩阵进行LU分解 /*for(i=0;i=500;i+) printf(%12le ,B0i);*/ for(k=1;k+)/开始迭代 a=0;for(i=0;i=500;i+) a=a+ui*ui;ita=sqrt(a);/计算2范数for(i=0;i=500;i+)/计算向量y(k-1) yi=ui/ita;for(i=0;i=500;i+) ui=0;d0=y0;/用LU分解解拟三角矩阵的方法求解u(k)for(i=1;it) t=i-1; else t=t; double temp1=0; for(;t=0;i-) if(i+2500) p=500; else p=i+2; double temp2=0; for(t=i+1;t=p;t+) temp2+=Bi-t+2t*ut; ui=(di-temp2)/B2i;for(i=0;i=500;i+) beta2=beta2+yi*ui;b=fabs(beta2-beta1);c=fabs(beta2);if(b/c)=sigma)/如果达到精度则结束并返回值 lam=1/beta2; goto finish;else beta1=beta2; beta2=0; finish:return(lam);/*矩阵分解程序*/void fenjie()/LU分解函数 int k,t,j,i,p; for(i=3;i=4;i+)/对于k=1时对应的矩阵元素先赋值 Bi0=Bi0/B20; for(k=2;k=501;k+)/从k=2开始将U中的元素存到系数矩阵中 if(k+2501) p=k+2;else p=501;for(j=k;jt) t=k-2; else t=t; if(j-2t) t=j-2; else t=t; double temp3=0; for(;t=k-1;t+) temp3+=Bk-t+2t-1*Bt-j+2j-1; Bk-j+2j-1=Bk-j+2j-1-temp3;if(k501) if(k+2501) p=k+2; else p=501; for(i=k+1;it) t=i-2; else t=t; if(k-2t) t=k-2; else t=t; double

温馨提示

  • 1. 本站所有资源如无特殊说明,都需要本地电脑安装OFFICE2007和PDF阅读器。图纸软件为CAD,CAXA,PROE,UG,SolidWorks等.压缩文件请下载最新的WinRAR软件解压。
  • 2. 本站的文档不包含任何第三方提供的附件图纸等,如果需要附件,请联系上传者。文件的所有权益归上传用户所有。
  • 3. 本站RAR压缩包中若带图纸,网页内容里面会有图纸预览,若没有图纸预览就没有图纸。
  • 4. 未经权益所有人同意不得将文件中的内容挪作商业或盈利用途。
  • 5. 人人文库网仅提供信息存储空间,仅对用户上传内容的表现方式做保护处理,对用户上传分享的文档内容本身不做任何修改或编辑,并不能对任何下载内容负责。
  • 6. 下载文件中如有侵权或不适当内容,请与我们联系,我们立即纠正。
  • 7. 本站不保证下载资源的准确性、安全性和完整性, 同时也不承担用户因使用这些下载资源对自己和他人造成任何形式的伤害或损失。

评论

0/150

提交评论