matlab程序总结.docx

1、matlab程序总结程序总结1、简单计算A:x255=xxxB:x255=xx2x4x8x16x32x64x1282. 算法1 (输入a(i)(i=0,1，,n),x;输出 算法2 （秦九韶算法）3.s=0;for i=1:n s=s+abs(x(i);End4.s=0; s=0;for i=1:n for i=1:n s=s+x(i)2; s=s+x(i)*x(i);end ends=sqrt(s) s=sqrt(s)5.s=0;for i=1:n if abs(x(i)s,s=abs(x(i);endend6. LU分解的matlap程序function A=lud(A)%功能：对方阵A作

2、三角分解A=LU,其中，% L为单位下三角阵，U为上三角阵，%输入：方阵A。%输出：紧凑存储A=LU.%注意：当A的主元=0时退出Matlabfor k=1:n-1for i=k+1:n if A(k,k) =0 quit; end A(i,k) =A(i,k)/ A(k,k); A(i,k+1:n)= A(i,k+1:n)- A(i,k) *A(k,k+1:n); endend7. 列主元Gauss消元法Lupd.m %功能：对方阵A作列主元三角分解PA=LU,其中，% L为单位下三角阵，U为上三角阵，排列阵P%用向量p表示。%输入：方阵A。%输出：紧凑存储LU=LU，以及p。%注意：当A奇

3、异时退出Matlab.function LU,p=lupd(A)%初始化n=length(A); p=1:n; LU=A;%分解过程for k=1:n%搜索列主元ik s,i=max(abs(LU(k:n,k); ik=i+k-1;%判断矩阵的奇异性if s=0 quit; end%行交换 if ik=k m=p(k); p(k)=p(ik); p(ik)=m; lk=LU(k,:); LU(k,:)=LU(ik,:); LU(ik,:)=lk; end%用消元法计算LU=LU if k=n break; end LU(k+1:n,k)=LU(k+1:n,k)/LU(k,k); LU(k+1:

4、n,k+1:n)=LU(k+1:n,k+1:n)-LU(k+1:n,k)* LU(k,k+1:n);End8. Householder矩阵变换function H,y=holder1(x)n=length(x);M=max(abs(x);if M=0, disp(M=0); return;else z=x/M;end;s=norm(z);if z(1)0 s=-s;endp=s*(s+z(1);u=z;u(1)= s+z(1);H=eye(n,n)-pu*u;y=zeros(n,1);y(1)=-M*s;9、解上三角方程function X=backsub(A,b)%InputA is an

5、nn upper- triangular nonsingullar matrix% -b is an n1 matrix%OutputX is the solution to the system AX=bn=length(b);X=zeros(n,1);X(n)=b(n)/A(n,n);for i=n-1:-1:1 X(i)=(b(i)-A(i,i+1:n)* X(i+1:n)/A(i,i);End10、matlap中高斯消元法function X=gauss(A,b)%InputA is an nn nonsingullar matrix% -b is an n1 matrix%Outpu

6、tX is the solution to the system AX=bn n=size(A); % 确定A的维数X=zeros(n,1);for k=1:n-1 for i=k+1:n % 消元过程 m=A(i,k)/ A(k,k); % A(k,k) 0 A(i,k+1:n)= A(i,k+1:n)-m*A(k,k+1:n); b(i)= b(i)-m*b(k); endendX=backsub(A, b); %回代求解11. 用矩阵分解法列主元的三角分解求解线性方程组lupdsv.m %功能：调用列主元三角分解函数LU,p=lupd(A)% 求解线性方程组Ax=b。%解法：PA=LU,

7、 Ax=bPAx=Pb % LUx=Pb, y=Ux % Ly=f=Pb, f(i)=b(p(i)%输入：方阵A，右端项b（行或列向量均可）%输出：解x（行向量）function x=lupdsv(A,b)n=length(b);LU,p=lupd(A);y(1)=b(p(1);for i=2:n y(i)=b(p(i)-LU(i,1:i-1)*y(1:i-1);endx(n)=y(n)/LU(n,n); for i=(n-1):-1:1 x(i)=(y(i)-LU(i,i+1:n)*x(i+1:n)/LU(i,i);end12. 用全主元的三角分解求解线性代数方程组function x=lu

8、pqdsv(A,b)n=length(b);LU,p,q=lupqd(A);y(1)=b(p(1);for i=2:n y(i)=b(p(i)-LU(i,1:i-1)*y(1:i-1);endz(n)=y(n)/LU(n,n);x(q(n)=z(n);for i=(n-1):-1:1 z(i)=(y(i)-LU(i,i+1:n)*z(i+1:n)/LU(i,i); x(q(i)=z(i);end13、G-S迭代法求解function x,k=gs(A,b)n n=size(A);x=zeros(1,n);for k=1:1000 error=0; for i=1:n s=0;xb=x(i);

9、for j=1:n if i=j,s=s+A(i,j)*x(j);end end x(i)=(b(i)-s)/A(i,i); error=error+abs(x(i)-xb); endif error/n0.0001,break;endendfprintf(k.no.=%3.0f,error=%7.2en,k,error)14.Gauss-Seidel迭代法参考程序: n=9;b(2:n+1,2:n+1)=0.02;U=zeros(n+2,n+2);e=0.000000001;for k=1:1000 %迭代求解 er=0; for j=2:n+1 for i=2:n+1 Ub=U(i,j);

10、 U(i,j)=(U(i-1,j)+U(i+1,j)+U(i,j-1)+U(i,j+1)+b(i,j)/4; er=er+abs(Ub-U(i,j); %估计当前误差end end if er/n2e,break;end %判断是否达到计算精度，如果达到则退出循环end15.追赶法程序function u,k,=kgs1(n)f=0.02*ones(n+2,n+2);a=-1*ones(1,n+2); b=4*ones(1,n+2);c=-1*ones(1,n+2); d=zeros(1,n+2);u=zeros(n+2,n+2);e=0.000000001;for k=1:1000 % 迭代

11、求解 er=0; for j=2:n+1 Ub=u(:,j); d(:)=f(:,j)+u(:,j-1)+u(:,j+1) ; % 块Gauss-Seidel迭代 z(2)=b(2);y(2,j)=d(2); for i=3:n+1 % 追赶法求解之追过程 求解Ly=d l(i)=a(i)/z(i-1); z(i)=b(i)-l(i)*c(i-1); y(i,j)=d(i)-l(i)*y(i-1,j); endu(n+1,j)=y(n+1,j)/z(n+1); % 追赶法求解之赶过程 求解Uz=y for m=n:-1:2 u(m,j)=(y(m,j)-c(m)*u(m+1,j)/z(m);

12、end er=er+norm(Ub-u(:,j),1); % 估计误差end if er/n2e ,break;end % 判断是否达到计算精度，如果达到则退出循环end function u,k=kgs2(n)f=2*1/(n+1)2*ones(n+2,n+2);a=-1*ones(1,n); b=4*ones(1,n);c=-1*ones(1,n); u=zeros(n+2,n+2);e=0.00001;for k=1:2000 er=0; for j=2:n+1 Ub=u(:,j); d(1:n)=f(2:n+1,j)+u(2:n+1,j-1)+u(2:n+1,j+1) ; x=zg(a

13、,b,c,d);u(2:n+1,j)=x; er=er+norm(Ub-u(:,j),1); end if er/n2e ,break;end end16.计算zclear;x=-8:0.5:8; y=x; X=ones(size(y)*x; Y=y*ones(size(x); R=sqrt(X.2+Y.2)+eps; Z=sin(R)./R; 17.计算方程的两个根函数定义行： function x1,x2=ff(a,b,c) H1 %ff.m:this file is to solve %algebra equation ax2+bx+c=0 帮助文本 %a,b,c:input argum

14、ents %x1,x2:output arguments,as are the %roots of equation dta=b2-4*a*c; % calculating % discriminant s1=-1*b+sqrt(dta); s2=-1*b-sqrt(dta);函数体 %calculation roots of equation x1=s1/2; x2=s2/2;18生成Hilbert矩阵 for i=1:3 for j=1:4 a(i,j)=1/(i+j-1); end end format rat19.计算矩阵中负数的个数i=0 ; c=0; x=-2 0 15 6 8 0 6 9 0 while ilength(x)+1 i=i+1 if x(i)0,c=c+1;end endc20用“求逆”法求解 方程的根tic xi=inv(A)*b; ti=toc eri=norm(x-xi) rei=norm(A*xi-b)/norm(b)（2）用“左除”法求解 tic;xd=Ab; td=toc, erd=norm(x-xd), red=norm(A*xd-b)/norm(b)