matlab使用法.docx
《matlab使用法.docx》由会员分享，可在线阅读，更多相关《matlab使用法.docx（22页珍藏版）》请在冰豆网上搜索。
matlab使用法
二分法
function[x0,k]=bisect1（fun1,a,b,ep）
ifnargin<4
ep=1e-5;
end
fa=feval（fun1,a）;
fb=feval（fun1,b）;
iffa*fb>0
x0=[fa,fb];
k=0;
return;
end
k=1;
whileabs（b-a）/2>ep
x=（a+b）/2;
fx=feval（fun1,x）;
iffx*fa<0
b=x;
fb=fx;
else
a=x;
fa=fx;
k=k+1;
end
end
x0=（a+b）/2;
>>fun1=inline（'x^3-x-1'）;
>>[x0,k]=bisect1（fun1,1.3,1.4,1e-4）
x0=
1.3247
k=
7
>>
简单迭代法
function[x0,k]=iterate1（fun1,x0,ep,N）
ifnargin<4
N=500;
end
ifnargin<3
ep=1e-5;
end
x=x0;
x0=x+2*ep;
k=0;
whileabs（x-x0）>ep&kx0=x;
x=feval（fun1,x0）;
k=k+1;
end
x0=x;
ifk==N
warning（'已达最大迭代次数'）
end
>>fun1=inline（'（x+1）^（1/3）'）;
>>[x0,k]=iterate1（fun1,1.5）
x0=
1.3247
k=
7
>>fun1=inline（'x^3-1'）;
>>[x0,k]=iterate1（fun1,1.5）
x0=
Inf
k=
9
>>
Steffesen加速迭代（简单迭代法的加速）
function[x0,k]=steffesen1（fun1,x0,ep,N）
ifnargin<4
N=500;
end
ifnargin<3
ep=1e-5;
end
x=x0;
x0=x+2*ep;
k=0;
whileabs（x-x0）>ep&kx0=x;
y=feval（fun1,x0）;
z=feval（fun1,y）;
x=x0-（y-x0）^2/（z-2*y+x0）;
k=k+1;
end
x0=x;
ifk==N
warning（'已达最大迭代次数'）
end
>>fun1=inline（'（x+1）^（1/3）'）;
>>[x0,k]=steffesen1（fun1,1.5）
x0=
1.3247
k=
3
>>fun1=inline（'x^3-1'）;
>>[x0,k]=steffesen1（fun1,1.5）
x0=
1.3247
k=
6
Newton迭代
function[x0,k]=Newton7（fname,dfname,x0,ep,N）
ifnargin<5
N=500;
end
ifnargin<4
ep=1e-5;
end
x=x0;
x0=x+2*ep;
k=0;
whileabs（x-x0）>ep&kx0=x;
x=x0-feval（fname,x0）/feval（dfname,x0）;
k=k+1;
end
x0=x;
ifk==N
warning（'已达最大迭代次数'）
end
>>fname=inline（'x-cos（x）'）;
>>dfname=inline（'1+sin（x）'）;
>>[x0,k]=Newton7（fname,dfname,pi/4,1e-8）
x0=
0.7391
k=4
非线性方程求根的Matlab函数调用举例：
1.求多项式的根：
求f（x）=x^3-x-1=0的根：
>>roots（[10-1-1]）
ans=
1.3247
-0.6624+0.5623i
-0.6624-0.5623i
2.求一般函数的根
>>fun=inline（'x*sin（x^2-x-1）','x'）
fun=
Inlinefunction:
fun（x）=x*sin（x^2-x-1）
>>fplot（fun,[-20.1]）;gridon
>>x=fzero（fun,[-2,-1]）
x=
-1.5956
>>x=fzero（fun,[-1-0.1]）
x=
-0.6180
[x,f,h]=fsolve（fun,-1.6）
x=
-1.5956
f=
1.4909e-009
h=
1
（h>0表示收敛，h<0表示发散，h=0表示已达到设定的计算函数值的最大次数）
第三章：
线性方程组的数值解法
1.高斯消元法
function[A,x]=gauss3（A,b）
%本算法用顺序高斯消元法求解线性方程组
n=length（b）;
A=[A,b];
fork=1:
n-1
A（（k+1）:
n,（k+1）:
（n+1））=A（（k+1）:
n,（k+1）:
（n+1））-A（（k+1）:
n,k）/A（k,k）*A（k,（k+1）:
（n+1））;
A（（k+1）:
n,k）=zeros（n-k,1）;
A;
end
x=zeros（n,1）;
%上面为消元过程
x（n）=A（n,n+1）/A（n,n）;
fork=n-1:
-1:
1
x（k）=（A（k,n+1）-A（k,（k+1）:
n）*x（（k+1:
n）））/A（k,k）;
end
%上面为回代过程
>>A=[234;352;4330];
>>b=[6,5,32]'
b=
6
5
32
>>[A,x]=gauss3（A,b）
A=
2.00003.00004.00006.0000
00.5000-4.0000-4.0000
00-2.0000-4.0000
x=
-13
8
2
列选主元的高斯消元法：
function[A,x]=gauss5（A,b）
%本算法用列选主元的高斯消元法求解线性方程组
n=length（b）;
A=[A,b];
fork=1:
n-1
%选主元
[ap,p]=max（abs（A（k:
n,k）））;
p=p+k-1;
ifp>k
t=A（k,:
）;
A（k,:
）=A（p,:
）;
A（p,:
）=t;
end
%消元
A（（k+1）:
n,（k+1）:
（n+1））=A（（k+1）:
n,（k+1）:
（n+1））-A（（k+1）:
n,k）/A（k,k）*A（k,（k+1）:
（n+1））;
A（（k+1）:
n,k）=zeros（n-k,1）;
end
%回代
x=zeros（n,1）;
x（n）=A（n,n+1）/A（n,n）;
fork=n-1:
-1:
1
x（k）=（A（k,n+1）-A（k,（k+1）:
n）*x（（sk+1:
n）））/A（k,k）;
end
>>A=[234;352;4330];b=[6,5,32]';
>>[A,x]=gauss5（A,b）
A=
4.00003.000030.000032.0000
02.7500-20.5000-19.0000
000.18180.3636
x=
-13
8
2
三角分解法：
Doolittle分解
function[L,U]=doolittle1（A）
n=length（A）;
U=zeros（n）;
L=eye（n）;
U（1,:
）=A（1,:
）;
L（2:
n,1）=A（2:
n,1）/U（1,1）;
fork=2:
n
U（k,k:
n）=A（k,k:
n）-L（k,1:
k-1）*U（1:
k-1,k:
n）;
L（k+1:
n,k）=A（k+1:
n,k）-L（k+1:
n,1:
k-1）*U（1:
k-1,n）/U（k,k）;
End
y=zeros（n,1）;
x=y;
y
（1）=b
（1）;
fori=2:
n
y（i）=b（i）-L（i,1:
i-1）*y（1:
i-1）;
end
x（n）=y（n）/U（n,n）;
fori=n-1:
-1:
1
x（i）=（y（i）-U（i,i+1:
n）*x（i+1:
n））/U（i,i）;
end
>>A=[123;252;315];b=[141820]';
>>[L,U,x]=doolittle1（A,b）
L=
100
210
3-81
U=
123
01-4
00-36
x=
2.8333
1.3333
2.8333
平方根法：
function[L,x]=choesky3（A,b）
n=length（A）;
L=zeros（n）;
L（:
1）=A（:
1）/sqrt（A（1,1））;
fork=2:
n
L（k,k）=A（k,k）-L（k,1:
k-1）*L（k,1:
k-1）';
L（k,k）=sqrt（L（k,k））;
fori=k+1:
n
L（i,k）=（A（i,k）-L（i,1:
k-1）*L（k,1:
k-1）'）/L（k,k）;
end
end
y=zeros（n,1）;
x=y;
y
（1）=b
（1）/L（1,1）;
fori=2:
n
y（i）=（b（i）-L（i,1:
i-1）*y（1:
i-1））/L（i,i）;
end
x（n）=y（n）/L（n,n）;
fori=n-1:
-1:
1
x（i）=（y（i）-L（i+1:
n,i）'*x（i+1:
n））/L（i,i）;
end
>>A=[4-11;-14.252.75;12.753.5]
A=
4.0000-1.00001.0000
-1.00004.25002.7500
1.00002.75003.5000
>>b=[467.25]'
b=
4.0000
6.0000
7.2500
[L,x]=choesky3（A,b）
L=
2.000000
-0.50002.00000
0.50001.50001.0000
x=
1
1
1
>>
迭代法求方程组的解
Jacobi迭代法：
function[x,k]=jacobi2（a,b,x0,ep,N）
%本算法用Jacobi迭代求解ax=b,用分量形式
n=length（b）;
k=0;
ifnargin<5
N=500;
end
ifnargin<4
ep=1e-5;
end
ifnargin<3
x0=zeros（n,1）;
y=zeros（n,1）;
end
x=x0;x0=x+2*ep;
whilenorm（x-x0,inf）>ep&kk=k+1;
x0=x;
fori=1:
n
y（i）=b（i）;
forj=1:
n
ifj~=i
y（i）=y（i）-a（i,j）*x0（j）;
end
end
ifabs（a（i,i））<1e-10|k==N
warning（'a（i,i）istoosmall'）;
return
end
y（i）=y（i）/a（i,i）;
end
x=y;
end
a=[430;34-1;0-14];b=[2430-24]'；
[x,k]=jacobi2（a,b）
x=
3.0000
4.0000
-5.0000
k=
59
Gauss-seidel迭代法：
function[x,k]=gaussseide2（a,b,x0,ep,N）
%本算法用Gauss-seidel迭代求解ax=b,用分量形式
n=length（b）;
k=0;
ifnargin<5
N=500;
end
ifnargin<4
ep=1e-5;
end
ifnargin<3
x0=zeros（n,1）;
y=zeros（n,1）;
end
x=x0;x0=x+2*ep;
whilenorm（x-x0,inf）>ep&kk=k+1;
x0=x;
y=x;
fori=1:
n
z（i）=b（i）;
forj=1:
n
ifj~=i
z（i）=z（i）-a（i,j）*x（j）;
end
end
ifabs（a（i,i））<1e-10|k==N
warning（'a（i,i）istoosmall'）;
return
end
z（i）=z（i）/a（i,i）;
x（i）=z（i）;
end
end
[x,k]=gaussseide2（a,b）
x=
3.0000
4.0000
-5.0000
k=
25
常微分方程数值解
function[x,y]=Euler1（fun,xspan,y0,h）
%本算法用欧拉格式计算微分方程y'=f（x,y）的解。
%fun为右端函数，xspan为求解的自变量区间，yo为初值，h为步长。
%输出向量x以及对应的函数值向量y.
x=xspan
（1）:
h:
xspan
（2）;
y
（1）=y0;
n=length（x）;
fork=1:
n-1
y（k+1）=y（k）+h*feval（fun,x（k）,y（k））;
end
x=x';y=y';
fun=inline（'y-2*x/y'）;
[x,y1]=Euler1（fun,[0,1],1,0.1）
x=
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
y1=
1
1.1
1.1918
1.2774
1.3582
1.4351
1.509
1.5803
1.6498
1.7178
1.7848
function[x,y]=Euler2（fun,xspan,y0,h）
%本算法用改进的欧拉格式计算微分方程y'=f（x,y）的解。
%fun为右端函数，xspan为求解的自变量区间，yo为初值，h为步长。
%输出向量x以及对应的函数值向量y.
x=xspan
（1）:
h:
xspan
（2）;
y
（1）=y0;
n=length（x）;
fork=1:
n-1
k1=feval（fun,x（k）,y（k））;
y（k+1）=y（k）+h*k1;
k2=feval（fun,x（k+1）,y（k+1））;
y（k+1）=y（k）+h*（k1+k2）/2;
end
x=x';y=y';
[x,y2]=Euler2（fun,[0,1],1,0.1）
x=
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
y2=
1
1.0959
1.1841
1.2662
1.3434
1.4164
1.486
1.5525
1.6165
1.6782
1.7379
y3=sqrt（1+2*x）
y3=
1
1.0954
1.1832
1.2649
1.3416
1.4142
1.4832
1.5492
1.6125
1.6733
1.7321
plot（x,y1,'+',x,y2,'o',x,y3,'r'）
将矩阵A按对角元绝对值由小到大顺序进行行列调动。
a=randn（4）；
a=
-0.4326-1.14650.3273-0.5883
-1.66561.19090.17462.1832
0.12531.1892-0.1867-0.1364
0.2877-0.03760.72580.1139
>>[u,v]=sort（-abs（diag（a）））
u=
-1.1909
-0.4326
-0.1867
-0.1139
v=
2
1
3
4
A=a（v,v）
A=
1.1909-1.66560.17462.1832
-1.1465-0.43260.3273-0.5883
1.18920.1253-0.1867-0.1364
-0.03760.28770.72580.1139
求二次多项式的根，避免上下溢，避免两个数相加或者相减接近0
function[xflag]=root2qe（a,b,c）
ifnargin==1
c=a（3）;b=a
（2）;a=a
（1）;
end
ifa~=0
M=max（abs（[a,b,c]））;
a=a/M;b=b/M;c=c/M;
s=sign（b）;
ifs==0
s=1;
end
x1=（-b-s*sqrt（b^2-4*a*c））/（2*a）;
x2=c/（a*x1）;
x=[x1x2];
flag='twosolutions';
elseifb~=0
x=-c/b;
flag='onesolution';
elseifc~=0
x=[];
flag='nosolution'
else
x=[1];
flag='allxaresolutions';
end
>>[xflag]=root2qe（6,10,-4）
x=
-2.00000.3333
flag=
twosolutions
>>[xflag]=root2qe（6e+154,10e+154,-4e+154）
x=
-2.00000.3333
flag=
twosolutions
>>[xflag]=root2qe（0,1,1）
x=
-1
flag=
onesolution
>>[xflag]=root2qe（1,-1e+5,1）
x=
1.0e+004*
10.00000.0000
flag=
twosolutions
计算如下的递推数列：
x1=1/3;x2=1/12;x（k+1）=2.25*x（k）-0.5*x（k-1）;可以证明，xk=4^（1-k）/3;是一个单调递减数列，但是计算数据如下：
functionex2
x
（1）=1/3;
x
（2）=1/12;
fork=3:
60
x（k）=2.25*x（k-1）-0.5*x（k-2）;
end
plot（x（1:
60））
想想是什么原因导致这种现象？
找出函数y=arctan（x）在区间【1,6】的11个等分点，做10次插值多项式。
打印出Newton形式的系数，在同一个图上作出在区间【0,8】的33个等分点上插值多项式的值与y=arctan（x）的值，由此得出什么结论？
functionex3
n=11;
x=linspace（1,6,n）;
h=（6-1）/（n-1）;
y=atan（x）';
forj=2:
n
y（1:
n+1-j,j）=diff（y（1:
n+2-j,j-1））/（（j-1）*h）;
end
y=y（1,:
）;
formatlonge
fprintf（'\nthecoefficientsofNewtoninterpolationis:
\n\n'）;
fprintf（'%12.9f',y）;
pz=[];
v=linspace（0,8,33）;
fort=v,
z=y（n）;
forj=n-1:
-1:
1
z=z*（t-x（j））+y（j）;
end
pz=[pzz];
end
plot（v,pz,'g*-',v,atan（v）,'b+:
'）;
fprintf（'\nxi,p（xi）,atan（xi）,error\n'）;
forj=1:
length（v）
fprintf（'%12.6f%12.6f%12.6f%12.6f',...
v（j）,pz（j）,atan（v（j））,abs（pz（j）-atan（v（j））））;
end
（注：
可编辑下载，若有不当之处，请指正，谢谢!
）
升级会员 