数值分析

第二章数值实验：运用追赶法求解三对角方程组的算法，求各段电路的电流量。1该电阻电路存在未知电流错误!未找到引用源。~错误!未找到引用源。，它们构成了特殊的8X8三对角矩阵。根据追赶法的计算原理，需编写一个用于计算追赶法的库函数：functionchase。2源程序functionx=Untitled2(a,b,c,f)%定义函数chasen=length(b);ifn-1==length(a)fori=n-1:-1:1a(i+1)=a(i);endend%将a设置为n维向量b(1)=b(1)-a(1)*c(1);c(1)=c(1)/b(1);f(1)=f(1)/b(1);fori=2:n-1c(i)=c(i)/b(i);b(i)=b(i)-a(i)*c(i-1);f(i)=(f(i)-a(i)*f(i-1))/b(i);endf(n)=(f(n)-a(n)*f(n-1))/(b(n)-a(n)*c(n-1));fori=(n-1):-1:1f(i)=f(i)-c(i)*f(i+1);endx=f;3创建关于电流i的8X8三对角矩阵A，并调用库函数chase（A,b），即可求出各个电流值。A=[2-2000000;-25-200000;0-25-20000;00-25-2000;000-25-200;0000-25-20;00000-25-2;000000-25];a=[-2-2-2-2-2-2-2];b=[25555555];c=[-2-2-2-2-2-2-2];f=[220/700000000];x=Untitled2(a,b,c,f)x=3.14271.57130.78550.39250.19570.09670.04600.0184A*x'ans=3.14290.0000-0.00000.00000.00000.0000-0.00004数据分析追赶法使计算公式简化，同时计算量减少，计算结果比较满意。但追赶法有使用的局限性，只对求解三对角方程组比较方便。第三章数值实验分别用Jacobi迭代法和Gauss-Seidel迭代法求解方程组，迭代初始向量取x0=[00000]。1、Jacobi迭代法源代码：%函数输入实参表中A为线性方程系数矩阵，b为常数向量，n为系数矩阵维数，x0为初始向量，e为自己定义的误差上线，N为允许的最大迭代次数，以防止死循环。functionjacobi(A,b,n,x0,e,N)I=eye(n);%构造生成迭代法中使用的B矩阵的逆矩阵以及g矩阵D=zeros(1,n);fori=1:nD(1,i)=A(i,i);endD1=1./D;D=diag(D);D1=diag(D1);B=I-D1*A;g=D1*b;count=0;error=1;whileerror=ex=B*x0+g;count=count+1;error=norm(x-x0);x0=x;if(count=N)disp('不收敛');return;endenddisp(x);计算结果A=[101234;19-12-3;2-173-5;32312-1;4-3-5-115];b=[12;-27;14;-17;12];n=5;e=1e-4;N=500;x0=[00000]';jacobi(A,b,n,x0,e,N)1.0002-2.00022.9997-1.99990.99982GAUSS-SEIDEL迭代法源程序function[x,n]=gauseidel(A,b,x0,exp,M)%求解线性方程组的迭代法，其中%A为方程组的系数矩阵%b为方程组的右端项%x0为迭代初始化向量%eps为精度要求，默认值为1e-5%M为最大迭代次数，默认值为100%x为方程组的解%n为迭代次数ifnargin==3exp=1.0e-6;M=200;elseifnargin==4M=200;elseifnargin==4M=200;elseifnargin3errorreturn;endD=diag(diag(A));%求A的对角矩阵L=-tril(A,-1);%求A的下三角阵U=-triu(A,1);%求A的上三角阵G=(D-L)\U;f=(D-L)\b;x=G*x0+f;n=1;%迭代次数whilenorm(x-x0)=expx0=x;x=G*x0+f;n=n+1;if(n=M)disp('迭代次数太多可能不收敛');return;endend运行结果A=[101234;1-9-12-3;2-173-5;32312-1;4-3-5-115];x0=[00000]';exp=1e-6;M=500;b=[12;-27;14;-17;12];[x,n]=gauseidel(A,b,x0,exp,M)x=1.0002-2.00022.9997-1.99990.9998n=164数据分析Gauss-Seidel迭代法比Jacobi迭代法效果好，收敛速度快Jacobi迭代法迭代次数增加迭代结果越接近精确结果，。第四章已知矩阵2912836-30112-147168-33636-42303663084-66190（1）试用幂法求按模最大的特征值与特征向量（2）用QR方法求矩阵的所有特征值幂法1源程序function[l,v,s]=powermethod(A,x0,eps)%幂法求模最大的特征值及对应的特征向量powermethod%输入:A:待求的矩阵%EPX:误差限%v0:初始向量%输出:l:求的矩阵的最大特征值%v:求的矩阵的主特征向量%s:迭代次数ifnargin==2eps=1e-6;endv=x0;M=5000;m=0;l=0;for(k=l:M)y=A*v;m=max(y);%m为按摩最大分量v=y/m;if(abs(m-l)eps)l=m;%到所需精度，退出，l为主特征值s=k;%s为迭代步数return;elseif(k==m)disp;l=m;s=M;elsel=M;endendend2计算结果A=[19066-8430;6630342-36;336-168147-112;30-3628291];x0=[1111]';[l,v,s]=powermethod(A,x0,eps)l=1.0000v=0.33331.0000-0.0000-0.5000s=344QR的方法源程序function[H,B]=Hessenberg(A)n=length(A);B=eye(n);fork=1:n-2;X=zeros(n-k,1)H=eye(n);fori=1:n-kX(i)=A(i+k,k);enda=max(abs(X));ifa==0.0breakendX=X/a;c=X(1);b1=sqrt(sum(X.^2));ifX(1)=0b1=-b1;endX(1)=X(1)-b1;b=b1^2-b1*c;H0=eye(n-k)-X*X'/b;fori=1:n-kforj=1:n-kH(i+k,j+k)=H0(i,j);endendA=H*A*H;B=B*H;end计算结果A=[19066-8430;6630342-36;336-168147-112;30-3628291];[H,B]=Hessenberg(A)H=1.000000001.00000000-0.93230.3618000.36180.9323B=1.00000000-0.19200.8797-0.43500-0.9775-0.21090.00510-0.08730.42620.9004所有特征值为1.000，-0.1920,-0.2109,0.9004数据分析幂法计算简便易行，但在特征值接近时或有N个线性无关的特征向量时，收敛速度慢，QR法目前是求解实矩阵的最有效的方法，适用于求所有矩阵的特征值和特征向量，但收敛的速度慢，需要改进采用原点位移法或Aiteken加快收敛速度。