matlab入门经典教程--第四章 数值计算

第四章数值计算

4.1引言

本章将花较大的篇幅讨论若干常见数值计算问题：线性分析、一元和多元函数分析、微积分、数据分析、以及常微分方程（初值和边值问题）求解等。但与一般数值计算教科书不同，本章的讨论重点是：如何利用现有的世界顶级数值计算资源MATLAB。至于数学描述，本章将遵循“最低限度自封闭”的原则处理，以最简明的方式阐述理论数学、数值数学和MATLAB计算指令之间的在联系及区别。

对于那些熟悉其他高级语言（如FORTRAN，Pascal，C++）的读者来说，通过本章，MATLAB 卓越的数组处理能力、浩瀚而灵活的M函数指令、丰富而友善的图形显示指令将使他们体验到解题视野的豁然开朗，感受到摆脱烦琐编程后的眉眼舒展。

对于那些经过大学基本数学教程的读者来说，通过本章，MATLAB精良完善的计算指令，自然易读的程序将使他们感悟“教程”数学的基础地位和局限性，看到从“理想化”简单算例通向科学研究和工程设计实际问题的一条途径。

对于那些熟悉MATLAB基本指令的读者来说，通过本章，围绕基本数值问题展开的容将使他们体会到各别指令的运用场合和在关系，获得综合运用不同指令解决具体问题的思路和借鉴。

由于MATLAB的基本运算单元是数组，所以本章容将从矩阵分析、线性代数的数值计算开始。然后再介绍函数零点、极值的求取，数值微积分，数理统计和分析，拟合和插值，Fourier分析，和一般常微分方程初值、边值问题。本章的最后讨论稀疏矩阵的处理，因为这只有在大型问题中，才须特别处理。

从总体上讲，本章各节之间没有依从关系，即读者没有必要从头到尾系统阅读本章容。读者完全可以根据需要阅读有关节次。除特别说明外，每节中的例题指令是独立完整的，因此读者可以很容易地在自己机器上实践。

MATLAB从5.3版升级到6.x版后，本章容的变化如下：

●MATLAB从6.0版起，其矩阵和特征值计算指令不再以LINPACK和EISPACK库为基础，

而建筑在计算速度更快、运行更可靠的LAPACK和ARPACK程序库的新基础上。因此，虽然各种矩阵计算指令没有变化，但计算结果却可能有某些不同。这尤其突出地表现在涉及矩阵分解、特征向量、奇异向量等的计算结果上。对此，用户不必诧异，因为构成空间的基向量时不唯一的，且新版的更可信。本书新版全部算例结果是在6.x版上给出的。

●在5.3版本中，泛函指令对被处理函数的调用是借助函数名字符串进行的。这种调用

方式在6.x版中已被宣布为“过渡期允许使用但即将被淘汰的调用方式”；而新的调用方式是借助“函数句柄”进行的。因此，关于述泛函指令，本章新版着重讲述如何使用“函数句柄”，同时兼顾“函数名字符串”调用法。

●MATLAB从6.0版起，提供了一组专门求微分方程“边值问题”数值解的指令。适应这

种变化，本章新增第4.14.5节，用2个算例阐述求解细节。

● 5.3版中的积分指令quad8已经废止；6.x版启用新积分指令quad l；6.5版新增三重

积分指令triplequad。本章新版对此作了相应的改变。

4.2LU分解和恰定方程组的解

4.2.1LU分解、行列式和逆

4.2.2恰定方程组的解

【例4.2.2-1】“求逆”法和“左除”法解恰定方程的性能对比

（1）

randn('state',0);

A=gallery('randsvd',100,2e13,2);

x=ones(100,1);

b=A*x;

cond(A)

ans =

1.9990e+013

（2）

tic

xi=inv(A)*b;

ti=toc

eri=norm(x-xi)

rei=norm(A*xi-b)/norm(b)

ti =

0.7700

eri =

0.0469

rei =

0.0047

（3）

tic;xd=A\b;

td=toc,erd=norm(x-xd),red=norm(A*xd-b)/norm(b)

td =

erd =

0.0078

red =

2.6829e-015

4.2.3数、条件数和方程解的精度

【例4.2.3-1】Hilbert矩阵是著名的病态矩阵。MATLAB中有专门的Hilbert矩阵及其准确

Hx 近似解和准确解进行比较。

逆矩阵的生成函数。本例将对方程b

N=[6 8 10 12 14];

for k=1:length(N)

n=N(k);

H=hilb(n);

Hi=invhilb(n);

b=ones(n,1);

x_approx=H\b;

x_exact=Hi*b;

ndb=norm(H*x_approx-b);nb=norm(b);

ndx=norm(x_approx - x_exact);nx=norm(x_approx);

er_actual(k)=ndx/nx;

K=cond(H);

er_approx(k)=K*eps;

er_max(k)=K*ndb/nb;

end

disp('Hilbert矩阵阶数'),disp(N)

format short e

disp('实际误差 er_actual'),disp(er_actual),disp('')

disp('近似的最大可能误差 er_approx'),disp(er_approx),disp('')

disp('最大可能误差 er_max'),disp(er_max),disp('')

Hilbert矩阵阶数

6 8 10 12 14

实际误差 er_actual

1.5410e-010 1.7310e-007 1.9489e-004 9.1251e-002

2.1257e+000

近似的最大可能误差 er_approx

3.3198e-009 3.3879e-006 3.5583e-003 3.9846e+000 9.0475e+001

最大可能误差 er_max

7.9498e-007 3.8709e-002 1.2703e+003 4.7791e+007 4.0622e+010

4.3矩阵特征值和矩阵函数

4.3.1特征值和特征向量的求取

【例4.3.1-1】简单实阵的特征值问题。

A=[1,-3;2,2/3];[V,D]=eig(A)

V =

0.7746 0.7746

0.0430 - 0.6310i 0.0430 + 0.6310i

D =

0.8333 + 2.4438i 0

0 0.8333 - 2.4438i

【例4.3.1-2】本例演示：如矩阵中有元素与截断误差相当时的特性值问题。A=[3 -2 -0.9 2*eps

-2 4 -1 -eps

-eps/4 eps/2 -1 0

-0.5 -0.5 0.1 1 ];

[V1,D1]=eig(A);ER1=A*V1-V1*D1

[V2,D2]=eig(A,'nobalance');ER2=A*V2-V2*D2

ER1 =

0.0000 0.0000 0.0000 0.0000

0 -0.0000 -0.0000 -0.0000

0.0000 -0.0000 -0.0000 0.0000

0.0000 0.0000 0.0000 -0.5216

ER2 =

1.0e-014 *

-0.2665 0.0111 -0.0559 -0.1055

0.4441 0.1221 0.0343 0.0833

0.0022 0.0002 0.0007 0

0.0194 -0.0222 0.0222 0.0333

【例4.3.1-3】指令eig与eigs的比较。

rand('state',1),A=rand(100,100)-0.5;

t0=clock;[V,D]=eig(A);T_full=etime(clock,t0)

options.tol=1e-8;

options.disp=0;

t0=clock;[v,d]=eigs(A,1,'lr',options);

T_part=etime(clock,t0)

[Dmr,k]=max(real(diag(D)));

d,D(1,1)

T_full =

0.2200

T_part =

3.1300

d =

3.0140 + 0.2555i

ans =

3.0140 + 0.2555i

vk1=V(:,k);

vk1=vk1/norm(vk1);v=v/norm(v);

V_err=acos(norm(vk1'*v))*180/pi

D_err=abs(D(k,k)-d)/abs(d)

V_err =

1.2074e-006

D_err =

4.2324e-010

4.3.2特征值问题的条件数

【例4.3.2-1】矩阵的代数方程条件数和特征值条件数。

B=eye(4,4);B(3,4)=1;B

format short e,c_equ=cond(B),c_eig=condeig(B)

B =

1 0 0 0

0 1 0 0

0 0 1 1

0 0 0 1

c_equ =

2.6180e+000

Warning: Matrix is close to singular or badly scaled.

Results may be inaccurate. RCOND = 1.110223e-016. > In D:\MATLAB6P1\toolbox\matlab\matfun\condeig.m at line 30 c_eig =

1.0000e+000

1.0000e+000

4.5036e+015

4.5036e+015

【例4.3.2-2】对亏损矩阵进行Jordan分解。

A=gallery(5)

[VJ,DJ]=jordan(A);

[V,D,c_eig]=condeig(A);c_equ=cond(A);

DJ,D,c_eig,c_equ

A =

-9 11 -21 63 -252

70 -69 141 -421 1684

-575 575 -1149 3451 -13801

3891 -3891 7782 -23345 93365

1024 -1024 2048 -6144 24572

DJ =

0 1 0 0 0

0 0 1 0 0

0 0 0 1 0

0 0 0 0 1

0 0 0 0 0

D =

Columns 1 through 4

-0.0408 0 0 0

0 -0.0119 + 0.0386i 0 0

0 0 -0.0119 - 0.0386i 0

0 0 0 0.0323 + 0.0230i

0 0 0 0

Column 5

0.0323 - 0.0230i

c_eig =

1.0e+010 *

2.1293

2.0796

2.0796

2.0020

2.0020

c_equ =

2.0253e+018

4.3.3复数特征值对角阵与实数块特征值对角阵的转化

【例4.3.3-1】把例4.3.1-1中的复数特征值对角阵D转换成实数块对角阵，使VR*DR/VR=A。A=[1,-3;2,2/3];[V,D]=eig(A);

[VR,DR]=cdf2rdf(V,D)

VR =

0.7746 0

0.0430 -0.6310

DR =

0.8333 2.4438

-2.4438 0.8333

4.3.4矩阵的谱分解和矩阵函数

【例4.3.4-1】数组乘方与矩阵乘方的比较。

clear,A=[1 2 3;4 5 6;7 8 9];

A_Ap=A.^0.3

A_Mp=A^0.3

A_Ap =

1.0000 1.2311 1.3904

1.5157 1.6207 1.7118

1.7928 1.8661 1.9332

A_Mp =

0.6962 + 0.6032i 0.4358 + 0.1636i 0.1755 - 0.2759i

0.6325 + 0.0666i 0.7309 + 0.0181i 0.8292 - 0.0305i

0.5688 - 0.4700i 1.0259 - 0.1275i 1.4830 + 0.2150i

【例4.3.4- 2】标量的数组乘方和矩阵乘方的比较。（A取自例4.3.4-1）pA_A=(0.3).^A

pA_M=(0.3)^A

pA_A =

0.3000 0.0900 0.0270

0.0081 0.0024 0.0007

0.0002 0.0001 0.0000

pA_M =

2.9342 0.4175 -1.0993

-0.0278 0.7495 -0.4731

-1.9898 -0.9184 1.1531

【例4.3.4-3】sin的数组运算和矩阵运算比较。（A取自例4.3.4-1）

A_sinA=sin(A)

A_sinM=funm(A,'sin')

A_sinA =

0.8415 0.9093 0.1411

-0.7568 -0.9589 -0.2794

0.6570 0.9894 0.4121

A_sinM =

-0.6928 -0.2306 0.2316

-0.1724 -0.1434 -0.1143

0.3479 -0.0561 -0.4602

4.4奇异值分解

4.4.1奇异值分解和矩阵结构

4.4.1.1奇异值分解定义

4.4.1.2矩阵结构的奇异值分解描述

4.4.2线性二乘问题的解

4.4.2.1矩阵除运算的广义化

4.4.2.2线性模型的最小二乘解

y ，进行三种解法比较。其中A取MATLAB库中的特【例4.4.2.2-1】对于超定方程Ax

殊函数生成。

（1）

A=gallery(5);A(:,1)=[];y=[1.7 7.5 6.3 0.83 -0.082]';

x=inv(A'*A)*A'*y,xx=pinv(A)*y,xxx=A\y

Warning: Matrix is close to singular or badly scaled.

Results may be inaccurate. RCOND = 5.405078e-018.

x =

3.4811e+000

5.1595e+000

9.5340e-001

-4.6569e-002

xx =

3.4759e+000

5.1948e+000

7.1207e-001

-1.1007e-001

Warning: Rank deficient, rank = 3 tol = 1.0829e-010.

xxx =

3.4605e+000

5.2987e+000

-2.9742e-001

（2）

nx=norm(x),nxx=norm(xx),nxxx=norm(xxx)

nx =

6.2968e+000

nxx =

6.2918e+000

nxxx =

6.3356e+000

（3）

e=norm(y-A*x),ee=norm(y-A*xx),eee=norm(y-A*xxx)

e =

6.9863e-001

ee =

4.7424e-002

eee =

4.7424e-002

4.5函数的数值导数和切平面

4.5.1法线

【例4.5.1-1】曲面法线演示。

y=-1:0.1:1;x=2*cos(asin(y));

[X,Y,Z]=cylinder(x,20);

surfnorm(X(:,11:21),Y(:,11:21),Z(:,11:21));

4.5.2偏导数和梯度

4.5.2.1理论定义

4.5.2.2数值计算指令

【例4.5.2.2-1】用一个简单矩阵表现diff和gradient指令计算方式。F=[1,2,3;4,5,6;7,8,9]

Dx=diff(F)

Dx_2=diff(F,1,2)

[FX,FY]=gradient(F)

[FX_2,FY_2]=gradient(F,0.5)

F =

1 2 3

4 5 6

7 8 9

Dx =

3 3 3

3 3 3

Dx_2 =

1 1

1 1

1 1

FX =

1 1 1

1 1 1

1 1 1

FY =

3 3 3

3 3 3

3 3 3

FX_2 =

2 2 2

2 2 2 2 2 2 FY_2 =

6 6 6 6 6 6 6 6 6

【例 4.5.2.2-2】研究偶极子(Dipole)的电势（Electric potential ）和电场强度（Electric field density ）。设在),(b a 处有电荷q +，在),(b a --处有电荷q -。那么在电荷所在平

面上任何一点的电势和场强分别为)1

1(4),(0-

+-=

r r q

y x V πε，V E -?= 。其中2222)()(,)()(b y a x r b y a x r +++=-+-=-+。

90

10941

?=πε。又设电荷6102-?=q ，5.1=a ，5.1-=b 。

clear;clf;q=2e-6;k=9e9;a=1.5;b=-1.5;x=-6:0.6:6;y=x; [X,Y]=meshgrid(x,y);

rp=sqrt((X-a).^2+(Y-b).^2);rm=sqrt((X+a).^2+(Y+b).^2); V=q*k*(1./rp-1./rm); [Ex,Ey]=gradient(-V);

AE=sqrt(Ex.^2+Ey.^2);Ex=Ex./AE;Ey=Ey./AE; cv=linspace(min(min(V)),max(max(V)),49); contourf(X,Y,V,cv,'k-') %axis('square')

title('\fontname{隶书}\fontsize{22}偶极子的场'),hold on quiver(X,Y,Ex,Ey,0.7) plot(a,b,'wo',a,b,'w+')

plot(-a,-b,'wo',-a,-b,'w-')

4.6 函数的零点

4.6.1多项式的根 4.6.2一元函数的零点

4.6.2.1

利用MATLAB 作图指令获取初步近似解 4.6.2.2 任意一元函数零点的精确解

【例 4.6.2.2-1】通过求t b e t t f at -=-)(sin )(2的零点，综合叙述相关指令的用法。 （1）

y=inline('sin(t)^2*exp(-a*t)-b*abs(t)','t','a','b');

%<1>

（2）

a=0.1;b=0.5;t=-10:0.01:10;

y_char=vectorize(y); % <3> Y=feval(y_char,t,a,b);

clf,plot(t,Y,'r');hold on,plot(t,zeros(size(t)),'k'); -10

-5

0510

-5-4

-3

-2

-1

1

t

y (t )

（3） 由于Notebook 中无法实现zoom 、ginput 指令涉及的图形和鼠标交互操作，因此下面指令必须在MATLAB 指令窗中运行，并得到如图4.6-2所示的局部放大图及鼠标操作线。

zoom on

[tt,yy]=ginput(5);zoom off

图 4.6-2

tt tt =

-2.0032 -0.5415 -0.0072

0.5876 1.6561

（4）

[t4,y4,exitflag]=fzero(y,tt(4),[],a,b)

%<11>

t4 =

0.5993 y4 = 0

exitflag = 1

（5）

[t3,y3,exitflag]=fzero(y,tt(3),[],a,b) t3 =

-0.5198 y3 =

-5.5511e-017 exitflag = 1

（6）

op=optimset('fzero') op =

ActiveConstrTol: [] ......

Display: 'notify' ......

TolX: 2.2204e-016 TypicalX: []

op=optimset('tolx',0.01); op.TolX ans = 0.0100

（7）

[t4n,y4n,exitflag]=fzero(y,tt(4),op,a,b) t4n = 0.6042 y4n = 0.0017 exitflag = 1

4.6.3多元函数的零点

【例 4.6.3-1】求解二元函数方程组???=+==-=0

)cos(),(0

)sin(),(21y x y x f y x y x f 的零点。

（1）

x=-2:0.5:2;y=x;[X,Y]=meshgrid(x,y); F1=sin(X-Y);F2=cos(X+Y); v=[-0.2, 0, 0.2]; contour(X,Y,F1,v)

（2）

[x0,y0]=ginput(2); disp([x0,y0])

-0.7926 -0.7843 0.7926 0.7843

（3）

fun='[sin(x(1)-x(2)),cos(x(1)+x(2))]'; %<12> [xy,f,exit]=fsolve(fun,[x0(2),y0(2)]) %<13> Optimization terminated successfully:

First-order optimality less than OPTIONS.TolFun, and no negative/zero curvature detected xy =

0.7854 0.7854 f =

1.0e-006 *

-0.0984 0.2011 exit = 1

〖说明〗

[fun.m]

function ff=fun(x) ff(1)=sin(x(1)-x(2)); ff(2)=cos(x(1)+x(2));

4.7 函数极值点

4.7.1一元函数的极小值点 4.7.2多元函数的极小值点

【例 4.7.2-1】求2

2

2)1()(100),(x x y y x f -+-=的极小值点。它即是著名的Rosenbrock's "Banana" 测试函数。该测试函数有一片浅谷，许多算法难以越过此谷。（演示本例搜索过程的文件名为exm04072_1_1.m 。） （1）

ff=inline('100*(x(2)-x(1)^2)^2+(1-x(1))^2','x');

（2）

x0=[-1.2,1];[sx,sfval,sexit,soutput]=fminsearch(ff,x0)

sx =

1.0000 1.0000

sfval =

8.1777e-010

sexit =

1

soutput =

iterations: 85

funcCount: 159

algorithm: 'Nelder-Mead simplex direct search'

（3）

[ux,sfval,uexit,uoutput,grid,hess]=fminunc(ff,x0)

Warning: Gradient must be provided for trust-region method;

using line-search method instead.

> In D:\MATLAB6P1\toolbox\optim\fminunc.m at line 211

Optimization terminated successfully:

Current search direction is a descent direction, and magnitude of

directional derivative in search direction less than 2*options.TolFun ux =

1.0000 1.0000

sfval =

1.9116e-011

uexit =

1

uoutput =

iterations: 26

funcCount: 162

stepsize: 1.2992

firstorderopt: 5.0020e-004

algorithm: 'medium-scale: Quasi-Newton line search'

grid =

1.0e-003 *

-0.5002

-0.1888

hess =

820.4028 -409.5496

-409.5496 204.7720

4.8数值积分

4.8.1一元函数的数值积分

4.8.1.1闭型数值积分

【例 4.8.1.1-1】求dx e I x ?

-=1

2

，其精确值为 74684204.0 。

（1）

syms x;IS=int('exp(-x*x)','x',0,1) vpa(IS) IS =

1/2*erf(1)*pi^(1/2) ans = .3185

（2）

fun=inline('exp(-x.*x)','x');

Isim=quad(fun,0,1),IL=quadl(fun,0,1) Isim =

0.7468 IL =

0.7468

（3）

Ig=gauss10(fun,0,1) Ig =

0.7463

（4）

xx=0:0.1:1.5;ff=exp(-xx.^2); pp=spline(xx,ff); int_pp=fnint(pp);

Ssp=ppval(int_pp,[0,1])*[-1;1] Ssp =

0.7468

（5）

图4.8-1

4.8.1.2 开型数值积分

[gauss10.m]

function g = gauss10(fun,a,b) %GAUSS10(fun,a,b) % fun

%====================================================== x = [0.1488743390;0.4333953941;0.6974095683;... 0.8650633667;0.9739065285];

w = [0.2955242247;0.2692667193;0.2190863625;... 0.1494513492;0.0666713443];

t = .5*(b+a)+.5*(b-a)*[-flipud(x);x]; W = [flipud(w);w];

g = sum(W.*feval(fun,t))*(b-a)/2;

【例 4.8.1.2-1】当)cos()(x x f =时，比较解析积分和近似积分。 （1）

syms x;F=int('cos(x)','x',-1,1),vpa(F) F =

2*sin(1) ans = 1.606 （2）

aF=cos(1/sqrt(3))+cos(-1/sqrt(3)) aF =

1.6758

【例4.8.1.2-2】求?

=1

1ln

dx x I ，准确结果是 .886226922

=π 。 （1）

syms x;IS=vpa(int('sqrt(log(1/x))','x',0,1))

Warning: Explicit integral could not be found.

> In D:\MATLAB6P5\toolbox\symbolic\sym\int.m at line 58 In D:\MATLAB6P5\toolbox\symbolic\char\int.m at line 9 IS = .7057

（2）用quad 指令求积

ff=inline('sqrt(log(1./x))','x');Isim=quad(ff,0,1) Warning: Divide by zero.

> In D:\MATLAB6P5\toolbox\matlab\funfun\inlineeval.m at line 13 In D:\MATLAB6P5\toolbox\matlab\funfun\inline\feval.m at line 34 In D:\MATLAB6P5\toolbox\matlab\funfun\quad.m at line 62 Isim =

0.8862

（3）

Ig=gauss10(ff,0,1) Ig =

0.8861

4.8.2多重数值积分

4.8.2.1 积分限为常数的二重积分指令

【例 4.8.2.1-1】计算?

???????=2

1

1001dy dx x S y x 和????

????=1

02112dy dx x S

y x 。

（1）

syms x y

ssx01=vpa(int(int(x^y,x,0,1),y,1,2)) ssx12=vpa(int(int(x^y,y,0,1),x,1,2))

Warning: Explicit integral could not be found.

> In D:\MATLAB6P5\toolbox\symbolic\sym\int.m at line 58 ssx01 =

ssx12 = 1.215

（2）

zz=inline('x.^y','x','y'); nsx01=dblquad(zz,0,1,1,2) nsx12=dblquad(zz,1,2,0,1) nsx01 = 0.4055 nsx12 =

1.2293

4.8.2.2 积分限为函数的二重积分

[double_int.m]

function SS=double_int(fun,innlow,innhi,outlow,outhi) %double_int % %fun %innlow,innhi %outlow,outhi

y1=outlow;y2=outhi;x1=innlow;x2=innhi;f_p=fun; SS=quad(G_yi,y1,y2,[],[],x1,x2,f_p);

[G_yi.m]

function f=G_yi(y,x1,x2,f_p) %G_yi %y %x1,x2 % %f_p

y=y(:);n=length(y);

if isnumeric(x1)==1;xx1=x1*ones(size(y));else xx1=feval(x1,y);end if isnumeric(x2)==1;xx2=x2*ones(size(y));else xx2=feval(x2,y);end for i=1:n

f(i)=quad(f_p,xx1(i),xx2(i),[],[],y(i)); end f=f(:);

【例 4.8.2.2-1】计算dy dx y x I y ?

???

????+=4

1

222)(。 （1）

[x_low.m]

function f=x_low(y) f=sqrt(y);

（2）

ff=inline('x.^2+y.^2','x','y'); SS=double_int(ff,x_low,2,1,4)

Warning: Minimum step size reached; singularity possible. > In D:\MATLAB6P5\toolbox\matlab\funfun\quad.m at line 88 In D:\MATLAB6p5\work\G_yi.m at line 11

In D:\MATLAB6P5\toolbox\matlab\funfun\inline\feval.m at line 20 In D:\MATLAB6P5\toolbox\matlab\funfun\quad.m at line 62 In D:\MATLAB6p5\work\double_int.m at line 8 SS =

9.5810

（4）

Ssym=vpa(int(int('x^2+y^2','x','sqrt(y)',2),'y',1,4)) Ssym = 9.524

4.8.3卷积

4.8.3.1 “完整”离散序列的数值卷积 4.8.3.2 “截尾”离散序列的数值卷积 4.8.3.3 多项式乘法与离散卷积的算法同构 4.8.3.4 连续时间函数的数值卷积 4.8.3.5 卷积的MATLAB 实现

【例4.8.3.5-1】有序列??

?==else

n n A 12

,,4,30

1

)( 和 ??

?==else

n n B 9

,,3,20

1

)( 。

（A ）

%

a=ones(1,10);n1=3;n2=12; b=ones(1,8);n3=2;n4=9;

c=conv(a,b);nc1=n1+n3;nc2=n2+n4; kc=nc1:nc2; %

aa=a(1:6);nn1=3;nn2=8;

cc=conv(aa,b);ncc1=nn1+n3; nx=nn2+n4;

ncc2=min(nn1+n4,nn2+n3);

kx=(ncc2+1):nx;kcc=ncc1:ncc2;N=length(kcc); stem(kcc,cc(1:N),'r','filled')

axis([nc1-2,nc2+2,0,10]),grid,hold on

stem(kc,c,'b'),stem(kx,cc(N+1:end),'g'),hold off

4

6

8

10

12

14

16

18

20

22

12345678910

【例4.8.3.5-2】求函数)()(t U e t u t

-=和)()(2

/t U te

t h t -=的卷积。本例展示：（A ）符号Laplace 变换求卷积的理论表示；（B ）SIMULINK 卷积法的执行过程和它的快速精确性。（C ）从理论符号解产生相应的理论数值序列。 （1）

syms tao;t=sym('t','positive'); US1=laplace(exp(-t)); HS1=laplace(t*exp(-t/2))

yt1=simple(ilaplace(US1*HS1)) HS1 =

1/(1/2+s)^2 yt1 =

4*exp(-t)+(2*t-4)*exp(-1/2*t)

（2）

图 4.8-3

（3）

t=yt2(:,1);

yyt1=eval(vectorize(char(yt1))); [dy,kd]=max(abs(yyt1-yt2(:,2))); dy12=dy/abs(yyt1(kd)) dy12 =

2.8978e-006

【例 4.8.3.5-3】用“零阶”近似法求)()(t U e t u t

-=和)()(2

/t U te t h t -=的卷积。本例演示：（A ）连续函数的有限长度采样。（B ）卷积数值计算三个误差（“截尾”误差、“零阶”近似误差、计算机字长误差）的影响。（C ）卷积“无截尾误差”区间、“非平凡”区间端点的确定。（D ）离散卷积和连续卷积之间的关系。（E ）指令conv 的使用。（F ）绘图分格线的运用。 （1）

（2）

%

t2=3;t4=11;T=0.01;

tu=0:T:t2;N2=length(tu);

th=0:T:t4;N4=length(th);

u=exp(-tu);h=th.*exp(-th/2);

tx=0:T:(t2+t4);Nx=length(tx);

yt3=T*conv(u,h);

%

t=tx;yyt1=eval(vectorize(char(yt1)));

[dy,kd]=max(abs(yyt1(1:N2)-yt3(1:N2)));

dy13(1)=dy/abs(yyt1(kd));

[dy,kd]=max(abs(yyt1(N2+1:N4)-yt3(N2+1:N4)));

dy13(2)=dy/abs(yyt1(N2+kd));

[dy,kd]=max(abs(yyt1(N4+1:Nx)-yt3(N4+1:Nx)));

dy13(3)=dy/abs(yyt1(N4+kd));

（3）

disp('与理论结果的相对误差')

disp([blanks(4),'[0,3]段 [3,11]段 [11,14]段']),disp(dy13) plot(t,yyt1,':b',tx,yt3,'r')

set(gca,'Xtick',[0,3,11,14]),grid

与理论结果的相对误差

[0,3]段 [3,11]段 [11,14]段

4.9随机数据的统计描述

4.9.1统计分布的数字特征

【例4.9.1-1】样本统计特征量计算示例。

%

X(:,1)=ones(10,1);X(1,1)=100;X(10,1)=0.01;

rand('state',1);randn('state',1)

X(:,2)=rand(10,1);

X(:,3)=randn(10,1);X(:,3)=2*abs(min(X(:,3)))+X(:,3);

%

Moment1.arithmetic=mean(X);Moment1.median=median(X); Moment1.geometric=geomean(X);Moment1.harmmonic=harmmean(X); %

Moment2.Standard=std(X);Moment2.variance=var(X);

Moment2.absolute=mad(X);Moment2.range=range(X);

%

X,Moment1,Moment2

X =

100.0000 0.9528 3.0699

1.0000 0.7041

2.2997

1.0000 0.9539 1.3535

1.0000 0.5982 3.0790

1.0000 0.8407 1.7674

1.0000 0.4428 1.7758

1.0000 0.8368 1.1027

1.0000 0.5187

2.6017

1.0000 0.0222 1.2405

0.0100 0.3759 2.3739

Moment1 =

arithmetic: [10.8010 0.6246 2.0664]

median: [1 0.6511 2.0377]

geometric: [1 0.4691 1.9463]

harmmonic: [0.0926 0.1682 1.8276]

Moment2 =

Standard: [31.3429 0.2951 0.7273]

variance: [982.3760 0.0871 0.5289]

absolute: [17.8398 0.2331 0.6184]

range: [99.9900 0.9317 1.9762]

4.9.2样本分布的频数直方图描述

【例4.9.2-1】hist指令的使用示例。

randn('state',1),rand('state',31)

x=randn(100,1);y=rand(100,1);

%

subplot(1,2,1),hist(x,7)

matlab入门经典范例

num1=[13]; den1=conv([1,1],[1,0]); G1=tf(num1,den1); num2=[5.096,13]; ssys1=conv([1,1],[1,0]); ssys2=conv([0.098,1],[1]); den2=conv( ssys1,ssys2); G2=tf(num2,den2); figure(1) margin(G1); hold on margin(G2); num1=[13]; den1=conv([1,1],[1,0]); G1=tf(num1,den1); num2=[5.096,13]; ssys1=conv([1,1],[1,0]); ssys2=conv([0.098,1],[1]); den2=conv( ssys1,ssys2); G2=tf(num2,den2); figure(1) margin(G1); hold on margin(G2); num=[4.56,10]; ssys1=conv([1,1],[1,0]); ssys2=conv([0.114,1],[1]); den=conv( ssys1,ssys2); G=tf(num,den); figure(1) bode(G) num=[4.56,10]; ssys1=conv([1,1],[1,0]); ssys2=conv([0.114,1],[1]); den=conv( ssys1,ssys2); G=tf(num,den); figure(1) msrgin(G); num1=[13]; den1=conv([1,1],[1,0]);

MATLAB经典数学建模教程

第 1 节Matlab 基本知识 一、Matlab 的主要功能 Matlab是一种功能非常强大的工程语言，诞生于20世纪70年代，1984年正式推向市场。2002年8月，Matlab6.5开始发布。是进行科学研究和产品开发必不可少的工具。 ●数值和符号计算 矩阵(数组)的四则运算（Matrix＋Laboratory）、数值差分、导数、积分、求解微分方程、微分方程的优化等 ●数字图像、数字信号处理 ●工程和科学绘图 ●控制系统设计 ●财务工程 ●建模、仿真功能 二、Matlab 的界面 1.命令窗口（Command Window）： Matlab各种操作命令都是由命令窗口开始，用户可以在命令窗口中输入Matlab命令，实现其相应的功能。此命令窗口主要包括文本的编辑区域和菜单栏（如：四则运算；“;”禁止显示变量的值；↑↓遍历以前的命令）。在命令窗口空白区域单击鼠标右键，打开快捷菜单，各项命令功能如下： Evaluate Selection :打开所选文本对应的表达式的值。 Open Selection :打开文本所对应的MatLab文件。 Cut :剪切编辑命令。 Paste :粘贴编辑命令。 2. M-文件编辑/调试（Editor/Debugger）窗口 Matlab Editor/Debugger窗口是一个集编辑与调试两种功能于一体的工具环境。 M-文件（函数文件） ●什么是M-文件：它是一种和Dos环境中的批处理文件相似的脚本文件，对于简单问题， 直接输入命令即可，但对于复杂的问题和需要反复使用的则需做成M－文件（Script File）。 ●创建M-文件的方法： Matlab命令窗的File/New/M-file。 在Matlab命令窗口运行edit。 ●M－文件的扩展名：*.m ●执行M-文件：F5 ●M文件的调试 选择Debug菜单，其各项命令功能如下： Step :逐步执行程序。 Step in :进入子程序中逐步执行调试程序。

matlab经典习题及解答

第1章 MATLAB概论 1.1与其他计算机语言相比较，MATLAB语言突出的特点是什么？ MATLAB具有功能强大、使用方便、输入简捷、库函数丰富、开放性强等特点。 1.2 MATLAB系统由那些部分组成？ MATLAB系统主要由开发环境、MATLAB数学函数库、MATLAB语言、图形功能和应用程序接口五个部分组成。 1.4 MATLAB操作桌面有几个窗口？如何使某个窗口脱离桌面成为独立窗口？又如何将脱离出去的窗口重新放置到桌面上？ 在MATLAB操作桌面上有五个窗口，在每个窗口的右上角有两个小按钮，一个是关闭窗口的Close按钮，一个是可以使窗口成为独立窗口的Undock按钮，点击Undock按钮就可以使该窗口脱离桌面成为独立窗口，在独立窗口的view菜单中选择Dock ……菜单项就可以将独立的窗口重新防止的桌面上。 1.5 如何启动M文件编辑/调试器？ 在操作桌面上选择“建立新文件”或“打开文件”操作时，M文件编辑/调试器将被启动。在命令窗口中键入edit命令时也可以启动M文件编辑/调试器。 1.6 存储在工作空间中的数组能编辑吗？如何操作？ 存储在工作空间的数组可以通过数组编辑器进行编辑：在工作空间浏览器中双击要编辑的数组名打开数组编辑器，再选中要修改的数据单元，输入修改内容即可。 1.7 命令历史窗口除了可以观察前面键入的命令外，还有什么用途？ 页脚内容1

命令历史窗口除了用于查询以前键入的命令外，还可以直接执行命令历史窗口中选定的内容、将选定的内容拷贝到剪贴板中、将选定内容直接拷贝到M文件中。 1.8 如何设置当前目录和搜索路径，在当前目录上的文件和在搜索路径上的文件有什么区别？ 当前目录可以在当前目录浏览器窗口左上方的输入栏中设置，搜索路径可以通过选择操作桌面的file 菜单中的Set Path菜单项来完成。在没有特别说明的情况下，只有当前目录和搜索路径上的函数和文件能够被MATLAB运行和调用，如果在当前目录上有与搜索路径上相同文件名的文件时则优先执行当前目录上的文件，如果没有特别说明，数据文件将存储在当前目录上。 1.9 在MATLAB中有几种获得帮助的途径？ 在MATLAB中有多种获得帮助的途径： （1）帮助浏览器：选择view菜单中的Help菜单项或选择Help菜单中的MATLAB Help菜单项可以打开帮助浏览器； （2）help命令：在命令窗口键入“help”命令可以列出帮助主题，键入“help 函数名”可以得到指定函数的在线帮助信息； （3）lookfor命令：在命令窗口键入“lookfor 关键词”可以搜索出一系列与给定关键词相关的命令和函数 （4）模糊查询：输入命令的前几个字母，然后按Tab键，就可以列出所有以这几个字母开始的命令和函数。 注意：lookfor和模糊查询查到的不是详细信息，通常还需要在确定了具体函数名称后用help命令显示详细信息。 第2章MATLAB矩阵运算基础 页脚内容2

BP神经网络matlab实例(简单而经典)

p=p1';t=t1'; [pn,minp,maxp,tn,mint,maxt]=premnmx(p,t); %原始数据归一化 net=newff(minmax(pn),[5,1],{'tansig','purelin'},'traingdx'); %设置网络,建立相应的BP网络 net.trainParam.show=2000; % 训练网络 net.trainParam.lr=0.01; net.trainParam.epochs=100000; net.trainParam.goal=1e-5; [net,tr]=train(net ,pn,tn); %调用TRAINGDM算法训练BP网络 pnew=pnew1'; pnewn=tramnmx(pnew,minp,maxp); anewn=sim(net,pnewn); %对BP网络进行仿真 anew=postmnmx(anewn,mint,maxt); %还原数据 y=anew'; 1、BP网络构建 （1）生成BP网络 = (,[1 2...],{ 1 2...},,,) net newff PR S S SNl TF TF TFNl BTF BLF PF R?维矩阵。 PR：由R维的输入样本最小最大值构成的2 S S SNl：各层的神经元个数。 [ 1 2...] TF TF TFNl：各层的神经元传递函数。 { 1 2...} BTF：训练用函数的名称。 （2）网络训练 = net tr Y E Pf Af train net P T Pi Ai VV TV [,,,,,] (,,,,,,) （3）网络仿真 = [,,,,] (,,,,) Y Pf Af E perf sim net P Pi Ai T {'tansig','purelin'},'trainrp' BP网络的训练函数 训练方法训练函数 梯度下降法traingd 有动量的梯度下降法traingdm 自适应lr梯度下降法traingda 自适应lr动量梯度下降法traingdx 弹性梯度下降法trainrp Fletcher-Reeves共轭梯度法traincgf

matlab入门教程

MATLAB入门教程 1．MATLAB的基本知识 1-1、基本运算与函数 在MATLAB下进行基本数学运算，只需将运算式直接打入提示号（>>）之後，并按入Enter键即可。例如： >> (5*2+1.3-0.8)*10/25 ans =4.2000 MATLAB会将运算结果直接存入一变数ans，代表MATLAB运算後的答案（Answer）并显示其数值於萤幕上。 小提示： ">>"是MATLAB的提示符号（Prompt），但在PC中文视窗系统下，由於编码方式不同，此提示符号常会消失不见，但这并不会影响到MATLAB的运算结果。 我们也可将上述运算式的结果设定给另一个变数x： x = (5*2+1.3-0.8)*10^2/25 x = 42 此时MATLAB会直接显示x的值。由上例可知，MATLAB认识所有一般常用到的加（+）、减（-）、乘（*）、除（/）的数学运算符号，以及幂次运算（^）。 小提示：MATLAB将所有变数均存成double的形式，所以不需经过变数宣告（Variable declaration）。MATLAB同时也会自动进行记忆体的使用和回收，而不必像C语言,必须由使用者一一指定.这些功能使的MATLAB易学易用，使用者可专心致力於撰写程式，而不必被软体枝节问题所干扰。 若不想让MATLAB每次都显示运算结果，只需在运算式最後加上分号（；）即可，如下例： y = sin(10)*exp(-0.3*4^2); 若要显示变数y的值，直接键入y即可： >>y y =-0.0045 在上例中，sin是正弦函数，exp是指数函数，这些都是MATLAB常用到的数学函数。 下表即为MATLAB常用的基本数学函数及三角函数： 小整理：MATLAB常用的基本数学函数 abs(x)：纯量的绝对值或向量的长度 angle(z)：复数z的相角(Phase angle)

matlab经典编程例题

以下各题均要求编程实现，并将程序贴在题目下方。 1．从键盘输入任意个正整数，以0结束，输出那些正整数中的素数。 clc;clear; zzs(1)=input('请输入正整数：');k=1; n=0;%素数个数 while zzs(k)~=0 flag=0;%是否是素数，是则为1 for yz=2:sqrt(zzs(k))%因子从2至此数平方根 if mod(zzs(k),yz)==0 flag=1;break;%非素数跳出循环 end end if flag==0&zzs(k)>1%忽略0和1的素数 n=n+1;sus(n)=zzs(k); end k=k+1; zzs(k)=input('请输入正整数：'); end disp(['你共输入了' num2str(k-1) '个正整数。它们是：']) disp(zzs(1:k-1))%不显示最后一个数0 if n==0 disp('这些数中没有素数！')%无素数时显示 else disp('其中的素数是：') disp(sus) end 2．若某数等于其所有因子（不含这个数本身）的和，则称其为完全数。编程求10000以内所有的完全数。 clc;clear;

wq=[];%完全数赋空数组 for ii=2:10000 yz=[];%ii的因子赋空数组 for jj=2:ii/2 %从2到ii/2考察是否为ii的因子 if mod(ii,jj)==0 yz=[yz jj];%因子数组扩展，加上jj end end if ii==sum(yz)+1 wq=[wq ii];%完全数数组扩展，加上ii end end disp(['10000以内的完全数为：' num2str(wq)])%输出 3．下列这组数据是美国1900—2000年人口的近似值（单位：百万）。 （1）若. 2c + = y+ 与试编写程序计算出上式中的a、b、c; 的经验公式为 t at bt y （2）若.bt 的经验公式为 y= 与试编写程序计算出上式中的a、b; y ae t （3）在一个坐标系下，画出数表中的散点图（红色五角星），c + =2中 ax bx y+拟合曲线图(蓝色实心线),以及.bt y=(黑色点划线)。 ae （4）图形标注要求:无网格线,横标注“时间t”，纵标注“人口数(百万)”，图形标题“美国1900—2000年的人口数据”。 （5）程序中要有注释，将你的程序和作好的图粘贴到这里。 clf;clc;clear %清除图形窗、屏幕、工作空间 t=1900:10:2000; y=[76 92 106 123 132 151 179 203 227 250 281]; p1=polyfit(t,y,2);%二次多项式拟合

matlab2007教程

第 1 章基础准备及入门本章有三个目的：一是讲述MATLAB正常运行所必须具备的基础条件；二是简明地介 绍MATLAB及其操作桌面Desktop的基本使用方法；三是全面介绍MATLAB的帮助系统。 本章的前两节讲述：MATLAB的正确安装方法和MATLAB 环境的启动。因为指令窗是MATLAB最重要的操作界面，所以本章用第 1.3、1.4 两节以最简单通俗的叙述、算例讲述指令窗的基本操作方法和规则。这部分内容几乎对MATLAB各种版本都适用。第1.5到第1.8节专门介绍MATLAB最常用的另五个交互界面：历史指令窗、当前目录浏览器、工作空间浏览器、数组编辑器、M文件编辑器。鉴于实际应用中，帮助信息和求助技能的重要性。本章专设第1.9节专门叙述MATLAB的帮助体系和求助方法。 作者建议：不管读者此前是否使用过MATLAB，都不要忽略本章。 1.1 MATLAB的安装和工具包选择 MATLAB只有在适当的外部环境中才能正常运行。因此，恰当地配置外部系统是保证MATLAB运行良好的先决条件。MATLAB本身可适应于许多机种和系统，如 PC机和 Unix 工作站等。但本节只针对我国使用最广的PC机系统给予介绍。 对PC机用户来说，常常需要自己安装MATLAB。MATLAB R2007a（即旧编号MATLAB7.4）版要求Win2000或WinXP平台。下面介绍从光盘上安装MATLAB的方法。一般说来，当MATLAB光盘插入光驱后，会自启动“安装向导”。假如自启动没有实现，那么可以在<我的电脑>或<资源管理器>中双击setup.exe应用程序，使“安装向导”启动。安装过程中出现的所有界面都是标准的，用户只要按照屏幕提示操作，如输入用户名、单位名、口令等就行。 在安装MATLAB.R2007a时，会出现一个界面，该界面上有两个选项：Typical和Custom。假如你不熟悉MATLAB，或假如你机器的硬盘的自由空间远大于3G，或假如你需要用到光盘上MATLAB的所有功能及工具包，那么你就点选“Typical”。否则，点选“Custom”。 在点选“Custom”后，会引出如图1.1-1的界面。你可以根据需要，在“Select products to install”栏中勾选相应的组件。注意：MATLAB软件光盘总包含很多工具包，它们有的是通用的，有的则专业性很强。对一般用户来说，完全不必采取全部安装，而应根据需要有所选择。否则将占据很多硬盘空间。表1.1-1对各组件的描述供用户选择时参考。

Matlab基础教程

1-1、基本运算与函数 在MATLAB下进行基本数学运算，只需将运算式直接打入提示号（>>）之後，并按入Enter键即可。例如： >> (5*2+1.3-0.8)*10/25 ans =4.2000 MATLAB会将运算结果直接存入一变数ans，代表MATLAB运算後的答案（Answer）并显示其数值於萤幕上。 小提示： ">>"是MATLAB的提示符号（Prompt），但在PC中文视窗系统下，由於编码方式不同，此提示符号常会消失不见，但这并不会影响到MATLAB的运算结果。 我们也可将上述运算式的结果设定给另一个变数x： x = (5*2+1.3-0.8)*10^2/25 x = 42 此时MATLAB会直接显示x的值。由上例可知，MATLAB认识所有一般常用到的加（+）、减（-）、乘（*）、除（/）的数学运算符号，以及幂次运算（^）。 小提示： MATLAB将所有变数均存成double的形式，所以不需经过变数宣告（Variable declaration）。MATLAB同时也会自动进行记忆体的使用和回收，而不必像C语言,必须由使用者一一指定.这些功能使的MATLAB易学易用，使用者可专心致力於撰写程式，而不必被软体枝节问题所干扰。 若不想让MATLAB每次都显示运算结果，只需在运算式最後加上分号（；）即可，如下例： y = sin(10)*exp(-0.3*4^2);

若要显示变数y的值，直接键入y即可： >>y y =-0.0045 在上例中，sin是正弦函数，exp是指数函数，这些都是MATLAB常用到的数学函数。 下表即为MATLAB常用的基本数学函数及三角函数： 小整理：MATLAB常用的基本数学函数 abs(x)：纯量的绝对值或向量的长度 angle(z)：复数z的相角(Phase angle) sqrt(x)：开平方 real(z)：复数z的实部 imag(z)：复数z的虚部 conj(z)：复数z的共轭复数 round(x)：四舍五入至最近整数 fix(x)：无论正负，舍去小数至最近整数 floor(x)：地板函数，即舍去正小数至最近整数 ceil(x)：天花板函数，即加入正小数至最近整数 rat(x)：将实数x化为分数表示 rats(x)：将实数x化为多项分数展开

matlab学习笔记(入门)

数据类：double，unit8，unit16，unit32，int8，int16，int32，single，char，logical！Matlab中所有数值计算都可以用double类来进行！，unit8实际中最常用的图像 图像类型：亮度图像，二值图像，索引图像，RGB图像 亮度图像：是数据矩阵，若是unit8或uint16则是【0,255】或者是【0，65535】，若是double 类，则像素取值是浮点数 二值图像只有：0和1的逻辑数组！ 、 简单操作： 读图并显示详细情况 >> f=imread('E:\image\book.pgm');whos Name Size Bytes Class Attributes f 289x338 97682 uint8 将图像垂直翻转： >> f=imread('E:\image\book.pgm');fp=f(end:-1:1, : );imshow(fp) 将图像上下左右翻转： f=imread('E:\image\book.pgm');fc=f(end:-1:1,end:-1:1);imshow(fc) 将图像二次采样并显示详情： >> fs=f(1:2:end,1:2:end);imshow(fs) >> whos fs Name Size Bytes Class Attributes fs 145x169 24505 uint8 将图像取出一部分： >> fg=f(200:250,200:300);imshow(fg) 显示图像中的一条水平扫描线： >> plot(f(200, : ) 将两幅图像进行相乘： f=imread('c:\image\liangdian.jpg');g=imread('c:\image\shuiguo.jpg'); g=g(300:715,500:1149);f=f(1:416,1:650);f d=double(f);gd=double(g); p=fd.*gd;数组乘！ pmax=max(p(:));pmin=min(p(:));取最大最小值！ pn=mat2gray(p);figure,imshow(pn) 亮度变化： 函数imadjust是对灰度图像进行亮度变化的基本ipt工具： g=imadjust（f，[low-in high-in]，[low-in high-in]，gamma） Gamma为1线性映射，大于1，则映射被加权至更低（更暗的）输出值，小于一，加权至更高的输出值 明暗反转图像（负片）参数不同： >> f=imread('E:\image\book.pgm');g=imadjust(f, [0 1],[1 0 ]);imshow(g) >> f=imread('E:\image\book.pgm');g=imadjust(f, [0 1],[1 0 ],2);imshow(g) >> f=imread('E:\image\book.pgm');g=imadjust(f, [0 1],[1 0 ],0.5);imshow(g) 另外也可以这样：进行明暗反转： g=imcomplement(f);imshow(g) 将0.5到0.75之间的灰度级拓展到0-1，可用于突出我们感兴趣的亮度带

matlab教程详解 (5)

第五章符号计算 符号计算的特点：一，运算以推理解析的方式进行，因此不受计算误差积累问题困扰；二，符号计算，或给出完全正确的封闭解，或给出任意精度的数值解（当封闭解不存在时）；三，符号计算指令的调用比较简单，经典教科书公式相近；四，计算所需时间较长，有时难以忍受。 在MATLAB中，符号计算虽以数值计算的补充身份出现，但涉及符号计算的指令使用、运算符操作、计算结果可视化、程序编制以及在线帮助系统都是十分完整、便捷的。 MATLAB的升级和符号计算内核Maple的升级，决定着符号计算工具包的升级。但从用户使用角度看，这些升级所引起的变化相当细微。即使这样，本章还是及时作了相应的更新和说明。如MATLAB 6.5+ 版开始启用Maple VIII的计算引擎，从而克服了Maple V计算“广义Fourier变换”时的错误（详见第5.4.1节）。 5.1符号对象和符号表达式 5.1.1符号对象的生成和使用 【例5.1.1-1】符号常数形成中的差异 a1=[1/3,pi/7,sqrt(5),pi+sqrt(5)] % <1> a2=sym([1/3,pi/7,sqrt(5),pi+sqrt(5)]) % <2> a3=sym([1/3,pi/7,sqrt(5),pi+sqrt(5)],'e') % <3> a4=sym('[1/3,pi/7,sqrt(5),pi+sqrt(5)]') % <4> a24=a2-a4 a1 = 0.3333 0.4488 2.2361 5.3777 a2 = [ 1/3, pi/7, sqrt(5), 6054707603575008*2^(-50)] a3 = [ 1/3-eps/12, pi/7-13*eps/165, sqrt(5)+137*eps/280, 6054707603575008*2^(-50)] a4 = [ 1/3, pi/7, sqrt(5), pi+sqrt(5)] a24 = [ 0, 0, 0, 189209612611719/35184372088832-pi-5^(1/2)] 【例5.1.1-2】演示：几种输入下产生矩阵的异同。 a1=sym([1/3,0.2+sqrt(2),pi]) % <1> a2=sym('[1/3,0.2+sqrt(2),pi]') % <2> a3=sym('[1/3 0.2+sqrt(2) pi]') % <3> a1_a2=a1-a2 % a1 = [ 1/3, 7269771597999872*2^(-52), pi] a2 = [ 1/3, 0.2+sqrt(2), pi] a3 = [ 1/3, 0.2+sqrt(2), pi] a1_a2 =

Matlab2012教程--经典教程

第1章基础准备及入门 1.1 最简单的计算器使用法 为易于学习，本节以算例方式叙述，并通过算例归纳一些MATLAB最基本的规则和语 法结构。建议读者，在深入学习之前，先读一读本节。 2 【例1.3-1】求[122(74)]3的算术运算结果。本例演示：最初步的指令输入形式 和必需的操作步骤。 （1）用键盘在MA TLAB指令窗中输入以下内容 >> (12+2*(7-4))/3^2 （2）在上述表达式输入完成后，按[Enter] 键，该指令被执行，并显示如下结果。 ans = 2 〖说明〗 本例在指令窗中实际运行的情况参见图 1.3-1。 指令行“头首”的“>>”是“指令输入提示符”，它是自动生成的。本书在此后的输入指令前将不再带提示符“>>”。理由是：（A）为使本书简洁；（B）本书用MATLAB 的M-book写成，而在M-book中运行的指令前是没有提示符的。 5

MATLAB的运算符（如+、- 等）都是各种计算程序中常见的习惯符号。 一条指令输入结束后，必须按[Enter] 键，那指令才被执行。 由于本例输入指令是“不含赋值号的表达式”，所以计算结果被赋给MATLAB的一个默认变量“ans”。它是英文“answer”的缩写。 【例1.3-2】“续行输入”法。本例演示：或由于指令太长，或出于某种需要，输入指令行必 须多行书写时，该如何处理。 S=1-1/2+1/3-1/4+ ... 1/5-1/6+1/7-1/8 S = 0.6345 〖说明〗 MA TLAB用3个或3个以上的连续黑点表示“续行”，即表示下一行是上一行的继续。 本例指令中包含“赋值号”，因此表达式的计算结果被赋给了变量S。 指令执行后，变量S被保存在MA TLAB 的工作空间（Workspace）中，以备后用。如果用户不用clear 指令清除它，或对它重新赋值，那么该变量会一直保存在工作空间中， 直到本MATLAB 指令窗被关闭为止。 1.3.3数值、变量和表达式 前节算例只是表演了“计算器”功能，那仅是MA TLAB全部功能中小小一角。为深入 学习MA TLAB，有必要系统介绍一些基本规定。本节先介绍关于变量的若干规定。 一数值的记述 MATLAB的数值采用习惯的十进制表示，可以带小数点或负号。以下记述都合法。 3 -99 0.001 9.456 1.3e-3 4.5e33 在采用IEEE浮点算法的计算机上，数值通常采用“占用64位内存的双精度”表示。 其相对精度是eps （MATLAB的一个预定义变量），大约保持有效数字16位。数值范围大308308 致从10到10。 二变量命名规则 变量名、函数名是对字母大小写敏感的。如变量myvar和MyVar表示两个不同的变量。 sin是MATLAB定义的正弦函数名，但SIN,Sin等都不是。 变量名的第一个字符必须是英文字母，最多可包含63个字符（英文、数字和下连符）。 如myvar201是合法的变量名。 变量名中不得包含空格、标点、运算符，但可以包含下连符。如变量名my_var_201是合法的，且读起来更方便。而my,var201由于逗号的分隔，表示的就不是一个变量名。 6

MATLAB课程设计报告(绝对完整)

课程设计任务书 学生姓名：董航专业班级：电信1006班 指导教师：阙大顺，李景松工作单位：信息工程学院 课程设计名称：Matlab应用课程设计 课程设计题目：Matlab运算与应用设计5 初始条件： 1.Matlab6.5以上版本软件； 2.课程设计辅导资料：“Matlab语言基础及使用入门”、“Matlab及在电子信息课程中的应 用”、线性代数及相关书籍等； 3.先修课程：高等数学、线性代数、电路、Matlab应用实践及信号处理类相关课程等。 要求完成的主要任务：（包括课程设计工作量及其技术要求，以及说明书撰写等具体要求） 1.课程设计内容：根据指导老师给定的7套题目，按规定选择其中1套完成； 2.本课程设计统一技术要求：研读辅导资料对应章节，对选定的设计题目进行理论分析， 针对具体设计部分的原理分析、建模、必要的推导和可行性分析，画出程序设计框图，编写程序代码（含注释），上机调试运行程序，记录实验结果（含计算结果和图表），并对实验结果进行分析和总结。具体设计要求包括： ①初步了解Matlab、熟悉Matlab界面、进行简单操作； ②MA TLAB的数值计算：创建矩阵矩阵运算、多项式运算、线性方程组、数值统计； ③基本绘图函数：plot, plot3, mesh, surf等，要求掌握以上绘图函数的用法、简单图形 标注、简单颜色设定等； ④使用文本编辑器编辑m文件，函数调用； ⑤能进行简单的信号处理Matlab编程； ⑥按要求参加课程设计实验演示和答辩等。 3.课程设计说明书按学校“课程设计工作规范”中的“统一书写格式”撰写，具体包括： ①目录； ②与设计题目相关的理论分析、归纳和总结； ③与设计内容相关的原理分析、建模、推导、可行性分析； ④程序设计框图、程序代码（含注释）、程序运行结果和图表、实验结果分析和总结； ⑤课程设计的心得体会（至少500字）； ⑥参考文献（不少于5篇）； ⑦其它必要内容等。 时间安排：1.5周（分散进行） 参考文献： [1](美)穆尔,高会生,刘童娜,李聪聪．MA TLAB实用教程(第二版) . 电子工业出版社,2010. [2]王正林,刘明．精通MA TLAB(升级版) .电子工业出版社,2011. [3]陈杰. MA TLAB宝典(第3版) . 电子工业出版社,2011. [4]刘保柱,苏彦华,张宏林. MA TLAB 7.0从入门到精通(修订版) . 人民邮电出版社,2010. 指导教师签名：年月日 系主任（或责任教师）签名：年月日

MATLAB基础及应用教程

第4章程序设计 在前面我们已经看到，MATLAB不但可以在命令窗直接输入命令并运行，而且还可以生成自己的程序文件，这就是我们通常说的一类以M为后缀的M文件，本章我们就来研究这类文件的形成方法。 M文件可分分为两大类，一是命令式M文件（也称为脚本文件，script），二是函数式M 文件（function）。两类文件的区别在于： （1）命令式文件可以直接运行，函数式文件不能直接运行，只能调用。 （2）命令式文件运行时没有输入输出参量，函数式文件在调用时需要进行输入输出参量设置。 （3）命令式文件运行中可以调用工作空间的数据，运行中产生的所有变量为全局变量。 （4）函数式文件不能调用工作空间的数据，运行中产生的所有变量为局部变量。命令式文件运行中产生的所有变量为全局变量，可以调用和存储到工作空间的数据。 4.1 MATLAB的程序文件-M文件 4.1.1 脚本文件（Scripts） 当我们需要在命令窗进行大量的命令集合运行时，直接从命令窗口输入比较麻烦，这时就可以将这些命令集合存放在一个脚本文件（Scripts）中，运行时只需要输入其文件名就可以自动执行这些命令集合。需要注意的是，脚本文件运行所产生的变量都驻留在MATLAB 的工作空间中，同时脚本文件也可以调用工作空间中的数据。因此，脚本文件所涉及的变量是全局变量。前几章所涉及到的M文件都是这类脚本文件。 编辑一个脚本文件可以直接在命令窗口的左上角打开编辑窗进行编辑。 4.1.2 函数文件（function） 函数式文件（function）的构成 （1）函数定义行： Function [输出参量]=gauss（输入参量） （2）： 完成函数的功能。 （3）函数说明。 （4）函数行注。 从上面构成的情况看，函数式文件实际上是完成输入参量与输出参量的转换，这样的转换是由函数文件名为gauss的文件来完成的。函数体的功能必须说明清楚输入参量与输出参量的关系。函数说明是用来解释该函数的功能的，函数行注是对程序行进行说明的。上面（1）和（2）是必须的。 【例4-1】分析下面函数文件。 %一个数列，任意项等于前两项之和，输入项数可以给出这个数列

北航 MATLAB教程答案(张志涌)

1 数字1.5e2，1.5e3 中的哪个与1500相同吗？ 1.5e3 2 请指出如下5个变量名中，哪些是合法的？ abcd-2 xyz_ 3 3chan a 变量 ABCDefgh 2、5是合法的。 3 在MATLAB 环境中，比1大的最小数是多少？ 1+eps 4 设 a = -8 , 运行以下三条指令，问运行结果相同吗？为什么？ w1=a^(2/3) w2=(a^2)^(1/3) w3=(a^(1/3))^2 w1 = -2.0000 + 3.4641i ;w2 = 4.0000 ;w3 =-2.0000 + 3.4641i 5 指令clear, clf, clc 各有什么用处？ clear 清除工作空间中所有的变量。 clf 清除当前图形。clc 清除命令窗口中所有显示。 第二章 1 说出以下四条指令产生的结果各属于哪种数据类型，是“双精度”对象，还是“符号”符号对象？ 3/7+0.1双; sym(3/7+0.1)符; sym('3/7+0.1') 符;; vpa(sym(3/7+0.1)) 符; 2 在不加专门指定的情况下，以下符号表达式中的哪一个变量被认为是自由符号变量. sym('sin(w*t)'),sym('a*exp(-X)'),sym('z*exp(j*th)') symvar(sym('sin(w*t)'),1) w a z 3 （1）试写出求三阶方程05.443 =-x 正实根的程序。注意：只要正实根，不要出现其他根。 （2）试求二阶方程022=+-a ax x 在0>a 时的根。 （1）reset(symengine) syms x positive solve(x^3-44.5) ans = (2^(2/3)*89^(1/3))/2 （2）求五阶方程02 2 =+-a ax x 的实根 syms a positive %注意：关于x 的假设没有去除 solve(x^2-a*x+a^2) Warning: Explicit solution could not be found. > In solve at 83 ans = [ empty sym ] syms x clear syms a positive solve(x^2-a*x+a^2) ans = a/2 + (3^(1/2)*a*i)/2 a/2 - (3^(1/2)*a*i)/2 4 观察一个数（在此用@记述）在以下四条不同指令作用下的异同。 a =@, b = sym( @ ), c = sym( @ ,' d ' ), d = sym( '@ ' ) 在此，@ 分别代表具体数值 7/3 , pi/3 , pi*3^(1/3) ；而异同通过vpa(abs(a-d)) , vpa(abs(b-d)) , vpa(abs(c-d))等来观察。 ● 理解准确符号数值的创建法。 ● 高精度误差的观察。 （1）x=7/3 x=7/3;a=x,b=sym(x),c=sym(x,'d'),d=sym('7/3'),

Matlab经典案例

1、三维曲线 >> t=0:pi/50:10*pi; >> plot3(sin(2*t),cos(2*t),t) >> axis square >> grid on 2、一窗口多图形 >> t=-2*pi:0.01:2*pi; >> subplot(3,2,1) >> plot(t,sin(t)) >> subplot(3,2,2) >> plot(t,cos(t)) >> subplot(3,2,3) >> plot(t,tan(t)) >> axis([-pi pi -100 100]) >> subplot(3,2,4) >> plot(t,cot(t)) >> axis([-pi pi -100 100]) >> subplot(3,2,5) >> plot(t,atan(t)) >> subplot(3,2,6) >> plot(t,acot(t)) 3、图形样式、标注、题字(也可以利用菜单直接 Insert) >> x=0:pi/20:2*pi; >> plot(x,sin(x),'b-.') >> hold on >> plot(x,cos(x),'r--') >> hold on >> plot(x,sin(x)-1,'g:')

>> hold on >> plot(x,cos(x)-1) >> xlabel('x'); >> xlabel('x轴'); >> ylabel('y轴'); >> title('图形样式、标注等'); >> text(pi,sin(pi),'x=\pi'); >> legend('sin(x)','cos(x)','sin(x)-1','cos(x)-1'); >> [x1,y1]=ginput(1) %利用鼠标定位查找线上某点的值x1 = 2.0893 y1 = -0.5000 >> gtext('x=2.5') %鼠标定位放置所需的值在线上 4、 >> fplot('[sin(x),cos(x),sqrt(x)-1]',[0 2*pi]) M文件：myfun.m 内容如下： function y=myfun(x) y(:,1)=sin(x); y(:,2)=cos(x); y(:,3)=x^(1/2)-1; 再运行：>> fplot('myfun',[0 2*pi]) 同样可以得到右图 5、 >> [x,y]=fplot('sin',[0 2*pi]); >> [x1,y1]=fplot('cos',[0 2*pi]); >> plot(x,y,'-r',x1,y1,'-.k') >> legend('y=sinx','y=cosx') 6、

(Matlab)SVM工具箱快速入手简易教程

SVM工具箱快速入手简易教程(by faruto) 一. matlab 自带的函数(matlab帮助文件里的例 子)[只有较新版本的matlab中有这两个SVM的函数] ===== svmtrain svmclassify =====简要语法规则==== svmtrain Train support vector machine classifier Syntax SVMStruct = svmtrain(Training, Group) SVMStruct = svmtrain(..., 'Kernel_Function', Kernel_FunctionValue, ...) SVMStruct = svmtrain(..., 'RBF_Sigma', RBFSigmaValue, ...) SVMStruct = svmtrain(..., 'Polyorder', PolyorderValue, ...) SVMStruct = svmtrain(..., 'Mlp_Params', Mlp_ParamsValue, ...) SVMStruct = svmtrain(..., 'Method', MethodValue, ...) SVMStruct = svmtrain(..., 'QuadProg_Opts', QuadProg_OptsValue, ...) SVMStruct = svmtrain(..., 'SMO_Opts', SMO_OptsValue, ...) SVMStruct = svmtrain(..., 'BoxConstraint', BoxConstraintValue, ...) SVMStruct = svmtrain(..., 'Autoscale', AutoscaleValue, ...) SVMStruct = svmtrain(..., 'Showplot', ShowplotValue, ...) --------------------- svmclassify Classify data using support vector machine Syntax Group = svmclassify(SVMStruct, Sample) Group = svmclassify(SVMStruct, Sample, 'Showplot', ShowplotValue) ============================实例研究==================== load fisheriris %载入matlab自带的数据[有关数据的信息可以自己到UCI查找,这是UCI的经典数据之一],得到的数据如下图:

matlab经典代码大全

哈哈哈 MATLAB 显示正炫余炫图：plot(x,y1,'* r',x,y2,'o b') 定义【0,2π】；t=0:pi/10:2*pi; 定义函数文件：function [返回变量列表]=函数名（输入变量列表） 顺序结构：选择结构 1）if-else-end语句 其格式为： if 逻辑表达式 程序模块1； else 程序模块2； End 图片读取：%选择图片路径 [filename, pathname] = ... uigetfile({'*.jpg';'*.bmp';'*.gif'},'选择图片'); %合成路径+文件名 str=[pathname,filename]; %为什么pathname和filename要前面出现的位置相反才能运行呢？？？%读取图片 im=imread(str); %使用图片 axes(handles.axes1); %显示图片 imshow(im); 边缘检测： global im str=get(hObject,'string'); axes (handles.axes1); switch str case ' 原图' imshow(im); case 'sobel' BW = edge(rgb2gray(im),'sobel'); imshow(BW); case 'prewitt' BW = edge(rgb2gray(im),'prewitt');

imshow(BW); case 'canny' BW = edge(rgb2gray(im),'canny'); imshow(BW);Canny算子边缘定位精确性和抗噪声能力效果较好，是一个折中方案 end; 开闭运算： se=[1,1,1;1,1,1;1,1,1;1,1,1]; %Structuring Element I=rgb2gray(im); imshow(I,[]);title('Original Image'); I=double(I); [im_height,im_width]=size(I); [se_height,se_width]=size(se); halfheight=floor(se_height/2); halfwidth=floor(se_width/2); [se_origin]=floor((size(se)+1)/2); image_dilation=padarray(I,se_origin,0,'both'); %Image to be used for dilation image_erosion=padarray(I,se_origin,256,'both'); %Image to be used for erosion %%%%%%%%%%%%%%%%%% %%% Dilation %%% %%%%%%%%%%%%%%%%%% for k=se_origin(1)+1:im_height+se_origin(1) for kk=se_origin(2)+1:im_width+se_origin(2) dilated_image(k-se_origin(1),kk-se_origin(2))=max(max(se+image_dilation(k-se_origin(1):k+halfh eight-1,kk-se_origin(2):kk+halfwidth-1))); end end figure;imshow(dilated_image,[]);title('Image after Dilation'); %%%%%%%%%%%%%%%%% %%% Erosion %%% %%%%%%%%%%%%%%%%% se=se'; for k=se_origin(2)+1:im_height+se_origin(2) for kk=se_origin(1)+1:im_width+se_origin(1) eroded_image(k-se_origin(2),kk-se_origin(1))=min(min(image_erosion(k-se_origin(2):k+halfwidth -1,kk-se_origin(1):kk+halfheight-1)-se)); end end figure;imshow(eroded_image,[]);title('Image after Erosion'); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%% Opening(Erosion first, then Dilation) %%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%