matlab程序总结.docx

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）