数值代数课设

《数值代数课设》由会员分享，可在线阅读，更多相关《数值代数课设（14页珍藏版）》请在装配图网上搜索。

1、1. 用matlab编写一个Gauss消去法的程序求解方程组= 程序： >> clear >> A=[2 1 0;3 4 -1;5 2 6] A = 2 1 0 3 4 -1 5 2 6 >> b=[1;0;1] b = 1 0 1 >> [L,U]=lu(A) L = 0.4000 0.0714 1.0000 0.6000 1.0000 0 1.0000 0 0

2、U = 5.0000 2.0000 6.0000 0 2.8000 -4.6000 0 0 -2.0714 >> c=L\b c = 1.0000 -0.6000 0.6429 >> x=U\c x = 0.8621 -0.7241 -0.3103 2、列主元Gauss消去法求解方程组= 程序： A=input(请输入线性方程组的增广矩阵A=); n=length(A)-1; x=zeros(n,1); aa=zeros(n,1)

3、; for j=1:n for i=1:(n+1) AA(j,i)=abs(A(j,i)); end end for k=1:(n-1) for i=k:n aa(i-(k-1))=AA(i,k); end for i=k:n if AA(i,k)==max(aa) break end end if AA(i,k)==0 break fprintf(方程组系数矩阵奇异\n); else for j=k:(n+1) jh=A(i,j); A(i,j)=A(k,j); A(k,j)=jh; end end fenzi=A(k,k); for j=

4、k:(n+1) A(k,j)=A(k,j)/fenzi; end for p=(k+1):n jj=A(p,k); for j=k:(n+1) A(p,j)=A(p,j)-jj*A(k,j); end end end if k==(n-1) x(n)=A(n,(n+1))/A(n,n); for i=(n-1):(-1):1 he=0; for j=(i+1):n he=he+A(i,j)*x(j); end x(i)=A(i,(n+1))-he; end end x 请输入线性方程组的增广矩阵A=[2 4 0 1;3 -5 1 1;1 7 8 1]

5、 x = 0.4032 0.0484 0.0323 3、用LU分解求解线形方程组= 程序：>> A=[2 1 -1 5;1 6 4 -1;3 2 7 5;0 -1 8 2] A = 2 1 -1 5 1 6 4 -1 3 2 7 5 0 -1 8 2 >> b=[8;4;2;1] b = 8 4 2 1 >> [L,U]=lu(A) L = 0.66

6、67 -0.0625 -0.6692 1.0000 0.3333 1.0000 0 0 1.0000 0 0 0 0 -0.1875 1.0000 0 U = 3.0000 2.0000 7.0000 5.0000 0 5.3333 1.6667 -2.6667 0 0 8.3125 1.5000

7、0 0 0 2.5038 >> [L,U,P]=lu(A) L = 1.0000 0 0 0 0.3333 1.0000 0 0 0 -0.1875 1.0000 0 0.6667 -0.0625 -0.6692 1.0000 U = 3.0000 2.0000 7.0000 5.0000 0 5.3333 1.666

8、7 -2.6667 0 0 8.3125 1.5000 0 0 0 2.5038 P = 0 0 1 0 0 1 0 0 0 0 0 1 1 0 0 0 >> y=L\(P*b) y = 2.0000 3.3333 1.6250 7.9624 >> x=U\y x = -5.3063

9、 2.3333 -0.3784 3.1802 4、用Cholesky分解求解线形方程组= 程序：>> clear >> A=[6 3 -2;3 5 1;-2 1 7] A = 6 3 -2 3 5 1 -2 1 7 >> b=[2;-1;1] b = 2 -1 1 >> L=chol(A) L = 2.4495 1.2247 -0.8165 0 1.8708 1.0690 0

10、 0 2.2783 >> Y=L\b Y = 0.8165 -1.0690 1.2332 >> X=L\Y X = 0.9541 -0.8807 0.5413 5、利用Jacobic迭代法求解方程组=，eps=1.010，= 程序： function [x, k, index]=Jacobi(A, b, eps, it_max) if nargin <4 it_max=100; end if nargin <3 eps; end n=length(A); k=0; x=zeros(n,1)

11、; y=zeros(n,1); index=1; while 1 for i=1:n y(i)=b(i); for j=1:n if j~=i y(i)=y(i)-A(i,j)*x(j); end end if abs(A(i,i))<1e-10 | k==it_max index=0; return; end y(i)=y(i)/A(i,i);

12、 end if norm(y-x,inf)> A=[0.98 -0.05 -0.02;-0.04 -0.9 0.07;-0.02 0.09 0.94]； >> b=[1;1;1]； >> [x,k,index]=Jacobi(A,b,1e-6,100)； x = 0.9904 -1.0628 1.1867 k = 6 index = 1 6、利用G-S迭代法求解方程组=，esp=1.010,=

13、 程序： function y=gaussseidel(A,b,x0,esp) D=diag(diag(A)); U=-triu(A,1); L=-tril(A,-1); G=(D-L)\U; f=(D-L)\b; y=G*x0+f; n=1; while norm(y-x0)>=esp&n<=1000 x0=y; y=G*x0+f; n=n+1; end y n答案： A=[10 3 1;2 -10 3;1 3 10]; b=[14;-5;14]; x0=[0;0;0]; esp=1.0e-3; gaussseidel(

14、A,b,x0,esp) y = 1.00003896860354 1.00002773079501 0.99998778390114 n = 6 ans = 1.00003896860354 1.00002773079501 0.99998778390114 7、利用SOR迭代法求解方程组 = Eps=1.0*10^-6,,取 程序： function [n,x]=sor(A,b,X,nm,w,ww) %用超松弛迭代法求解方程组Ax=b %输入：A为方程组的系数矩阵，b为方程组右端

15、的列向量，X为迭代初值构成的列向量，nm为最大迭代次数，w为误差精度，ww为松弛因子 %输出：x为求得的方程组的解构成的列向量，n为迭代次数 n=1; m=length(A); D=diag(diag(A)); %令A=D-L-U,计算矩阵D L=tril(-A)+D; %令A=D-L-U,计算矩阵L U=triu(-A)+D; %令A=D-L-U,计算矩阵U M=inv(D-ww*L)*((1-ww)*D+ww*U); %计算迭代矩阵 g=ww*inv(D-ww*L)*b; %计算迭代格式

16、中的常数项 %下面是迭代过程 while n<=nm x=M*X+g; %用迭代格式进行迭代 if norm(x-X,inf)

17、： >> A=[0.78 -0.02 -0.12 -0.14;-0.02 0.86 -0.04 0.06;-0.12 -0.04 0.72 -0.08;-0.14 0.06 -0.08 0.74]; b=[0.76;0.08;1.12;0.68]; c=1000; d=1e-6; f=1.03; k=[0;0;0;0]; g=sor(A,b,k,c,d,f) 迭代次数为 n = 7 方程组的解为 x = 1.5350 0.1220 1.9752 1.4130 g =

18、7 8、用共轭梯度法求解上题的线性方程组， Eps=1.0*10^-6, 9、利用幂法求：的按模最大的特征值和对应的特征向量，限定最大的迭代步骤n=500，Eps=1.0，= 程序： function [k,lambda,Vk,Wc]=mifa(A,x0,eps,n) lambda=0;k=1;Wc =1; state=1; V=x0; while((k<=n)&(state==1)) Vk=A*V; [m j]=max(abs(Vk)); mk=m; tzw=abs(lambda-mk); Vk=(1/mk)*Vk; Txw=norm(V-Vk); Wc=ma

19、x(Txw,tzw); V=Vk;lambda=mk;state=0; if(Wc>eps) state=1; end k=k+1;Wc=Wc; end if(Wc<=eps) disp(迭代次数k,主特征值的近似值lambda,主特征向量的近似向量Vk,相邻两次迭代的误差Wc如下：) else disp(迭代次数k已经达到最大迭代次数max1,主特征值的迭代值lambda,主特征向量的迭代向量Vk,相邻两次迭代的误差Wc如下：) end Vk=V;k=k-1;Wc; 主程序： >> A=[1 2 3;2 3 4;3 4 5]; x0=[1;1;1]

20、; esp=1.0e-6; n=500; [k,lambda,xk,Wc]=mifa(A,x0,eps,n), [k,lambda,xk,Wc]=mifa(A,x0,0.000001,500), [x,D] = eig (A), Dzd=max(diag(D)), wuD= abs(Dzd- lambda), wux=x(:,2)./xk, 答案： 迭代次数k已经达到最大迭代次数max1,主特征值的迭代值lambda,主特征向量的迭代向量Vk,相邻两次迭代的误差Wc如下： k = 500 lambda = 9.6235 xk =

21、 0.5247 0.7623 1.0000 Wc = 1.7764e-015 迭代次数k,主特征值的近似值lambda,主特征向量的近似向量Vk,相邻两次迭代的误差Wc如下： k = 8 lambda = 9.6235 xk = 0.5247 0.7623 1.0000 Wc = 1.5190e-007 x = 0.8277 0.4082 0.3851 0.1424 -0.81

22、65 0.5595 -0.5428 0.4082 0.7339 D = -0.6235 0 0 0 -0.0000 0 0 0 9.6235 Dzd = 9.6235 wuD = 9.2421e-009 wux = 0.7781 -1.0710 0.4082 10、利用反幂法求A=的按模最小的特征值和对应的特征向量，限定最大的迭代步骤n=5

23、00，Eps=1.0，= 程序： function [k,lambdan,Vk,Wc]=ydwyfmf(A,V0,jlamb,jd,max1) [n,n]=size(A); A1=A-jlamb*eye(n); jd= jd*0.1;RA1=det(A1); if RA1==0 disp(请注意：因为A-aE的n阶行列式hl等于零，所以A-aE不能进行LU分解.) return end lambda=0; if RA1~=0 for p=1:n h(p)=det(A1(1:p, 1:p)); end hl=h(1:n); for i=1:n i

24、f h(1,i)==0 disp(请注意：因为A-aE的r阶主子式等于零，所以A-aE不能进行LU分解.) return end end if h(1,i)~=0 disp(请注意：因为A-aE的各阶主子式都不等于零，所以A-aE能进行LU分解.) k=1;Wc =1;state=1; Vk=V0; while((k<=max1)&(state==1)) [L U]=lu(A1); Yk=L\Vk;Vk=U\Yk; [m j]=max(abs(Vk)); mk=m;Vk1=Vk/mk; Yk1=L\Vk1;Vk1=U\Yk1; [m j]=max(

25、abs(Vk1)); mk1=m;Vk2=(1/mk1)*Vk1;tzw1=abs((mk-mk1)/mk1); tzw2=abs(mk1-mk);Txw1=norm(Vk)-norm(Vk1); Txw2=(norm(Vk)-norm(Vk1))/norm(Vk1); Txw=min(Txw1,Txw2); tzw=min(tzw1,tzw2); Vk=Vk2; mk=mk1; Wc=max(Txw,tzw); Vk=Vk2;mk=mk1;state=0; if(Wc>jd) state=1; end k=k+1;%Vk=Vk2,mk=mk1, end if(Wc<=jd

26、) disp(A-aE的秩R(A-aE)和各阶顺序主子式值hl、迭代次数k,按模最小特征值的近似值lambda,特征向量的近似向量Vk,相邻两次迭代的误差Wc如下：) else disp(A-aE的秩R(A-aE)和各阶顺序主子式值hl、迭代次数k已经达到最大迭代次数max1,按模最小特征值的迭代值lambda,特征向量的迭代向量Vk,相邻两次迭代的误差Wc如下：) end hl,RA1 end end [V,D]=eig(A,nobalance),Vk;k=k-1;Wc;lambdan=jlamb+1/mk1; 主程序： A=[1 2 0;4 -3 2;5

27、1 4]; x0=[1;1;1]; [k,lambda,Vk,Wc]=ydwyfmf(A,x0,0.2,0.000001,500) 请注意：因为A-aE的各阶主子式都不等于零，所以A-aE能进行LU分解. A-aE的秩R(A-aE)和各阶顺序主子式值hl、迭代次数k,按模最小特征值的近似值lambda,特征向量的近似向量Vk,相邻两次迭代的误差Wc如下： hl = 0.8000 -10.5600 -21.7280 RA1 = -21.7280 V = -0.1577 -0.5531 0.3776 -0.3243 -0.0508 -1.0000 -1.0000 1.0000 -0.1070 D = 5.1127 0 0 0 1.1836 0 0 0 -4.2964 k = 7 lambda = 1.1836 Vk = 0.5531 0.0508 -1.0000 Wc = 5.3678e-009