matlab程序总结.docx
- 文档编号:3583327
- 上传时间:2022-11-24
- 格式:DOCX
- 页数:10
- 大小:81.97KB
matlab程序总结.docx
《matlab程序总结.docx》由会员分享,可在线阅读,更多相关《matlab程序总结.docx(10页珍藏版)》请在冰豆网上搜索。
matlab程序总结
程序总结
1、简单计算A:
x255=x·x···x
B:
x255=x·x2·x4·x8·x16·x32·x64·x128
2.
算法1(输入a(i)(i=0,1,…,n),x;输出算法2(秦九韶算法)
3.
s=0;
fori=1:
n
s=s+abs(x(i));
End
4.
s=0;s=0;
fori=1:
nfori=1:
n
s=s+x(i)^2;s=s+x(i)*x(i);
endend
s=sqrt(s)s=sqrt(s)
5.
s=0;
fori=1:
n
ifabs(x(i))>s,s=abs(x(i));
end
end
6.LU分解的matlap程序
functionA=lud(A)
%功能:
对方阵A作三角分解A=LU,其中,
%L为单位下三角阵,U为上三角阵,
%输入:
方阵A。
%输出:
紧凑存储A=[L\U].
%注意:
当A的主元=0时退出Matlab
fork=1:
n-1
fori=k+1:
n
ifA(k,k)==0quit;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);
end
end
7.列主元Gauss消元法
Lupd.m
%功能:
对方阵A作列主元三角分解PA=LU,其中,
%L为单位下三角阵,U为上三角阵,排列阵P
%用向量p表示。
%输入:
方阵A。
%输出:
紧凑存储LU=[L\U],以及p。
%注意:
当A奇异时退出Matlab.
function[LU,p]=lupd(A)
%初始化
n=length(A);p=1:
n;LU=A;
%分解过程
fork=1:
n
%搜索列主元ik
[s,i]=max(abs(LU(k:
n,k)));ik=i+k-1;
%判断矩阵的奇异性
ifs==0quit;end
%行交换
ifik~=k
m=p(k);p(k)=p(ik);p(ik)=m;
lk=LU(k,:
);LU(k,:
)=LU(ik,:
);LU(ik,:
)=lk;
end
%用消元法计算LU=[L\U]
ifk==nbreak;end
LU(k+1:
n,k)=LU(k+1:
n,k)/LU(k,k);
LU(k+1:
n,k+1:
n)=LU(k+1:
n,k+1:
n)-LU(k+1:
n,k)*…
LU(k,k+1:
n);
End
8.Householder矩阵变换
function[H,y]=holder1(x)
n=length(x);
M=max(abs(x));
ifM==0,
disp(‘M=0');
return;
else
z=x/M;
end;
s=norm(z);
ifz
(1)<0
s=-s;
end
p=s*(s+z
(1));
u=z;
u
(1)=s+z
(1);
H=eye(n,n)-p\u*u';
y=zeros(n,1);
y
(1)=-M*s;
9、解上三角方程
functionX=backsub(A,b)
%Input—Aisann×nupper-triangularnonsingullarmatrix
%---bisann×1matrix
%Output—XisthesolutiontothesystemAX=b
n=length(b);
X=zeros(n,1);
X(n)=b(n)/A(n,n);
fori=n-1:
-1:
1
X(i)=(b(i)-A(i,i+1:
n)*X(i+1:
n))/A(i,i);
End
10、matlap中高斯消元法
functionX=gauss(A,b)
%Input—Aisann×nnonsingullarmatrix
%---bisann×1matrix
%Output—XisthesolutiontothesystemAX=b
[nn]=size(A);%确定A的维数
X=zeros(n,1);
fork=1:
n-1
fori=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);
end
end
X=backsub(A,b);%回代求解
11.用矩阵分解法列主元的三角分解求解线性方程组
lupdsv.m
%功能:
调用列主元三角分解函数[LU,p]=lupd(A)
%求解线性方程组Ax=b。
%解法:
PA=LU,Ax=b←→PAx=Pb
%LUx=Pb,y=Ux
%Ly=f=Pb,f(i)=b(p(i))
%输入:
方阵A,右端项b(行或列向量均可)
%输出:
解x(行向量)
functionx=lupdsv(A,b)
n=length(b);
[LU,p]=lupd(A);
y
(1)=b(p
(1));
fori=2:
n
y(i)=b(p(i))-LU(i,1:
i-1)*y(1:
i-1)';
end
x(n)=y(n)/LU(n,n);
fori=(n-1):
-1:
1
x(i)=(y(i)-LU(i,i+1:
n)*x(i+1:
n)')/LU(i,i);
end
12.用全主元的三角分解求解线性代数方程组
functionx=lupqdsv(A,b)
n=length(b);
[LU,p,q]=lupqd(A);
y
(1)=b(p
(1));
fori=2:
n
y(i)=b(p(i))-LU(i,1:
i-1)*y(1:
i-1)';
end
z(n)=y(n)/LU(n,n);x(q(n))=z(n);
fori=(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);
end
13、G-S迭代法求解
function[x,k]=gs(A,b)
[nn]=size(A);
x=zeros(1,n);
fork=1:
1000
error=0;
fori=1:
n
s=0;xb=x(i);
forj=1:
n
ifi~=j,s=s+A(i,j)*x(j);end
end
x(i)=(b(i)-s)/A(i,i);
error=error+abs(x(i)-xb);
end
iferror/n<0.0001,break;end
end
fprintf('k.no.=%3.0f,error=%7.2e\n',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;
fork=1:
1000%迭代求解
er=0;
forj=2:
n+1
fori=2:
n+1
Ub=U(i,j);
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
ifer/n^2 end 15.追赶法程序 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; fork=1: 1000%迭代求解 er=0; forj=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); fori=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); end u(n+1,j)=y(n+1,j)/z(n+1);%追赶法求解之赶过程求解Uz=y form=n: -1: 2 u(m,j)=(y(m,j)-c(m)*u(m+1,j))/z(m); end er=er+norm(Ub-u(: j),1);%估计误差 end ifer/n^2 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; fork=1: 2000 er=0; forj=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,b,c,d);u(2: n+1,j)=x'; er=er+norm(Ub-u(: j),1); end ifer/n^2 end 16.计算z clear;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: thisfileistosolve %algebraequationax^2+bx+c=0 帮助文本%a,b,c: inputarguments %x1,x2: outputarguments,asarethe %rootsofequation dta=b^2-4*a*c;%calculating %discriminant s1=-1*b+sqrt(dta); s2=-1*b-sqrt(dta); 函数体%calculationrootsofequation x1=s1/2; x2=s2/2; 18生成Hilbert矩阵 fori=1: 3 forj=1: 4 a(i,j)=1/(i+j-1); end end formatrat 19.计算矩阵中负数的个数 i=0;c=0;x=[-20156–80690] whilei i=i+1 ifx(i)<0,c=c+1;end end c 20用“求逆”法求解方程的根 tic xi=inv(A)*b; ti=toc eri=norm(x-xi) rei=norm(A*xi-b)/norm(b) (2)用“左除”法求解 tic;xd=A\b; td=toc, erd=norm(x-xd), red=norm(A*xd-b)/norm(b)
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- matlab 程序 总结