matlab程序总结Word格式文档下载.docx
《matlab程序总结Word格式文档下载.docx》由会员分享，可在线阅读，更多相关《matlab程序总结Word格式文档下载.docx（10页珍藏版）》请在冰豆网上搜索。
5.
ifabs（x（i））>
s,s=abs（x（i））;
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
7.列主元Gauss消元法
Lupd.m
对方阵A作列主元三角分解PA=LU,其中，
%L为单位下三角阵，U为上三角阵，排列阵P
%用向量p表示。
紧凑存储LU=[L\U]，以及p。
当A奇异时退出Matlab.
function[LU,p]=lupd（A）
%初始化
n=length（A）;
p=1:
n;
LU=A;
%分解过程
%搜索列主元ik
[s,i]=max（abs（LU（k:
n,k）））;
ik=i+k-1;
%判断矩阵的奇异性
ifs==0quit;
%行交换
ifik~=k
m=p（k）;
p（k）=p（ik）;
p（ik）=m;
lk=LU（k,:
）;
LU（k,:
）=LU（ik,:
LU（ik,:
）=lk;
%用消元法计算LU=[L\U]
ifk==nbreak;
LU（k+1:
n,k）=LU（k+1:
n,k）/LU（k,k）;
n,k+1:
n）=LU（k+1:
n）-LU（k+1:
n,k）*…
LU（k,k+1:
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）<
s=-s;
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）;
10、matlap中高斯消元法
functionX=gauss（A,b）
nnonsingullarmatrix
[nn]=size（A）;
%确定A的维数
fori=k+1:
n%消元过程
m=A（i,k）/A（k,k）;
%A（k,k）≠0
n）-m*A（k,k+1:
b（i）=b（i）-m*b（k）;
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）
[LU,p]=lupd（A）;
y
（1）=b（p
（1））;
fori=2:
y（i）=b（p（i））-LU（i,1:
i-1）*y（1:
i-1）'
x（n）=y（n）/LU（n,n）;
fori=（n-1）:
x（i）=（y（i）-LU（i,i+1:
n）*x（i+1:
n）'
）/LU（i,i）;
12.用全主元的三角分解求解线性代数方程组
functionx=lupqdsv（A,b）
[LU,p,q]=lupqd（A）;
z（n）=y（n）/LU（n,n）;
x（q（n））=z（n）;
z（i）=（y（i）-LU（i,i+1:
n）*z（i+1:
x（q（i））=z（i）;
13、G-S迭代法求解
function[x,k]=gs（A,b）
x=zeros（1,n）;
1000
error=0;
fori=1:
xb=x（i）;
forj=1:
ifi~=j,s=s+A（i,j）*x（j）;
x（i）=（b（i）-s）/A（i,i）;
error=error+abs（x（i）-xb）;
iferror/n<
0.0001,break;
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;
1000%迭代求解
er=0;
forj=2:
n+1
fori=2:
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
ifer/n^2<
e,break;
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）;
1000%迭代求解
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）;
u（n+1,j）=y（n+1,j）/z（n+1）;
%追赶法求解之赶过程求解Uz=y
form=n:
2
u（m,j）=（y（m,j）-c（m）*u（m+1,j））/z（m）;
er=er+norm（Ub-u（:
j）,1）;
%估计误差
e,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;
2000
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'
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矩阵
3
4
a（i,j）=1/（i+j-1）;
formatrat
19.计算矩阵中负数的个数
i=0;
c=0;
x=[-20156–80690]
whilei<
length（x）+1
i=i+1
ifx（i）<
0,c=c+1;
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）
升级会员 