MATLAB具有以下几个特点: 易学、适用范围广、功能强、开放性强、网络资源丰富。 启动 点击MATLAB图标,进入到MATLAB命令窗(Matlab Command Window)。 学会使用help命令。 学会使用demo命令。说明其功能强大。 演示 census; spinner;truss; pend.m plot([-0.2,0.2],[0;0],'color','y','linestyle','-','linewidth',10); g=0.98;l=1; theta0=pi/6;x0=l*sin(theta0); y0=-l*cos(theta0); axis([-0.75,0.75,-1.25,0]); axis('off'); head=line(x0,y0,'color','r','linestyle','.','erasemode','xor','markersize',40); body=line([0;x0],[0,y0],'color','b','linestyle','-','erasemode','xor'); t=0; dt=0.01; while t<=50 t=t+dt; theta=theta0*cos(sqrt(g/l)*t); x=l*sin(theta);y=-l*cos(theta); set(head,'xdata',x,'ydata',y); set(body,'xdata',[0;x],'ydata',[0;y]); drawnow; end 退出 在工具栏中点击File按钮,在下拉式菜单中单击Exit MATLAB项即可。 或者,在指令窗内键入exit或quit亦可。 矩阵运算的操作(demo) MATLAB的符号运算功能 求和 symsum(S) 对通项S求和,其中k为变量。且从0变到k-1。 symsum(S,v) 对通项S求和,指定其中v为变量。且v从0变到v-1。 symsum(S,a,b) 对通项S求和,其中k为变量。且从a变到b。 symsum(S,v,a,b) 对通项S求和,指定其中v为变量。且v从a变到b。 例:键入k=sym('k');symsum(k) 得 ans = 1/2*k^2-1/2*k 又例如:键入 symsum(k^2,0,10)得 ans = 385 又例如:键入symsum('x'^k/sym('k!'),k,0,inf)得 ans = exp(x) 这最后的一个例子是无穷项求和。 ⅱ 求导数 diff(S,v) 求表达式S对变量v的一阶导数。 diff(S,v,n) 求表达式S对变量v的n阶导数。 例如:键入命令 A=sym('[1/(1+a),(b+x)/cos(x);1,exp(x^2)]'); diff(A,'x') 得 ans = [ 0, 1/cos(x)+(b+x)/cos(x)^2*sin(x)] [ 0, 2*x*exp(x^2)] 又如求sin(x)+ex的三阶导数,键入命令 diff('sin(x)+x*exp(x)',3) 得 ans = -cos(x)+3*exp(x)+x*exp(x) 再如:求 A = [ x*sin(y), x^n+y] [ 1/x/y, exp(i*x*y)] 的先对x再对y的混合偏导数。 可键入命令: S=sym('[x*sin(y),x^n+y;1/x/y,exp(i*x*y)]'); dsdxdy=diff(diff(S,'x'),'y') 得: dsdxdy = [ cos(y), 0] [ 1/x^2/y^2, i*exp(i*x*y)-y*x*exp(i*x*y)] 求y=(lnx)x的导数 可键入命令: p='(log(x))^x'; p1=diff(p,'x') 得 p1 = log(x)^x*(log(log(x))+1/log(x)) 求y=xf(x2)的导数 可键入命令: p='x*f(x^2)'; p1=diff(p,'x') 得 p1 = f(x^2)+2*x^2*D(f)(x^2) 求xy=ex+y的导数 可键入命令: p='x*y(x)-exp(x+y(x))'; p1=diff(p,'x') p1 = y(x)+x*diff(y(x),x)-(1+diff(y(x),x))*exp(x+y(x)) p2='y+x*dy-(1+dy)*exp(x+y)=0'; dy=solve(p2,'dy')%把dy作为变量解方程 得 dy= -(y-exp(x+y))/(x-exp(x+y)) ⅲ 求极限 limit(P) 表达式P中自变量趋于零时的极限。 limit(P,a) 表达式P中自变量趋于a时的极限。 limit(P,x,a,'left') 表达式P中自变量x趋于a时的左极限。 limit(P,x,a,'right')表达式P中自变量x趋于a时的右极限。 例如:键入 P=sym('sin(x)/x'); limit(P) 得 ans = 1 键入 P=sym('1/x'); limit(P,x,0,'right') 得 ans = inf 键入 P=sym('(sin(x+h)-sin(x))/h');h=sym('h'); limit(P,h,0) 得 ans = cos(x) 键入 v=sym('[(1+a/x)^x,exp(-x)]'); limit(v,x,inf,'left') 得 ans = [ exp(a), 0] ⅳ 求泰勒展开式 taylor(f,v) f对v的五阶Maclaurin展开。 taylor(f,v,n) f对v的n-1阶Maclaurin展开。 例如求sin(x)e-x 的7阶Maclaurin展开。可键入 f=sym('sin(x)*exp(-x)');F=taylor(f,8) 得 F = x-x^2+1/3*x^3-1/30*x^5+1/90*x^6-1/630*x^7 如果要求sin(x)e-x 在x=1 处的7阶Taylor展开。可键入 f=sym('sin(x)*exp(-x)');F=taylor(f,8,1) 得 F = sin(1)*exp(-1)+(-sin(1)*exp(-1)+cos(1)*exp(-1))*(x-1) -cos(1)*exp(-1)*(x-1)^2 +(1/3*sin(1)*exp(-1)+1/3*cos(1)*exp(-1))*(x-1)^3 -1/6*sin(1)*exp(-1)*(x-1)^4 +(1/30*sin(1)*exp(-1)-1/30*cos(1)*exp(-1))*(x-1)^5 +1/90*cos(1)*exp(-1)*(x-1)^6 +(-1/630*cos(1)*exp(-1)-1/630*sin(1)*exp(-1))*(x-1)^7 多元函数的taylor展开 MATLAB不能直接进行多元函数的taylor展开。必须先调用MAPLE函数库中的mtaylor命令。方法为: 在MATLAB的工作窗口中键入 maple('readlib(mtaylor)') mtaylor的格式为 mtaylor(f,v,n) f为欲展开的函数式。 v为变量名。写成向量的形式:[var1=p1,var2=p2,…,varn=pn],展开式将在(p1,p2,…,pn)处进行。如只有变量名,将在0点处展开。n为展开式的阶数(n-1阶)。要完成taylor展开,只需键入maple('mtaylor(f,v,n)')即可。 例:在(x0,y0,z0)处将F=sin(x,y,z)进行2阶taylor展开。键入 syms x0 y0 z0 maple('readlib(mtaylor)'); maple('mtaylor(sin(x*y*z),[x=x0,y=y0,z=z0],2)') 得: ans = sin(x0*y0*z0)+cos(x0*y0*z0)*y0*z0*(x-x0)+cos(x0*y0*z0)*x0*z0*(y-y0)+cos(x0*y0*z0)*x0*y0*(z-z0) ⅴ 求积分 int(P) 对表达式P进行不定积分。 int(P,v) 以v为积分变量对P进行不定积分。 int(P,v,a,b) 以v为积分变量,以a为下限,b为上限对P进行定积分。 例如可键入int('-2*x/(1+x^2)^2') 得 ans = 1/(1+x^2) 键入int('x/(1+z^2)','z') 得 ans = atan(z)*x 键入int('x*log(1+x)',0,1) 得 ans = 1/4 定积分的上下限可以是(符号)函数。例如可键入: int('2*x','sin(t)','log(t)') 得 ans = log(t)^2-sin(t)^2 对(符号)矩阵进行积分,例 输入int('[exp(t),exp(a*t)]'),得: ans = [ exp(t), 1/a*exp(a*t)] ⑶ 求符号方程的解 ⅰ线性方程组的求解 线性方程组的形式为A*X=B;其中A至少行满秩。 X=linsolve(A,B) 输出方程的特解X。 例如:键入 A=sym('[cos(t),sin(t);sin(t),cos(t)]'); B=sym('[1;1]'); c=linsolve(A,B) c = [ 1/(sin(t)+cos(t))] [ 1/(sin(t)+cos(t))] 例如:键入 a=sym('[2,7,3,1;3,5,2,2;9,4,1,7]');b=sym('[6;4;2]'); X=linsolve(a,b) Warning: System is rank deficient. Solution is not unique. X = [ 0] [ 0] [ 2] [ 0] ⅱ 代数方程的求解 solve(P,v) 对方程P中的指定变量v求解。v可省略。 solve(p1,P2,…,Pn,v1,v2,…,vn) 对方程P1,P2,…Pn中的指定变量v1, v2…vn求解。 例:可输入 solve('p+sin(x)=r') 得: ans = -asin(p-r) 又例:可输入: P1='x^2+x*y+y=3';P2='x^2-4*x+3=0'; [x,y]=solve(P1,P2) 得: x = [ 1] [ 3] y = [ 1] [ -3/2] 可输入: P1='a+u^2+v^2=0';P2='u-v=1'; [u,v]=solve(P1,P2,'u','v') 得: u = [ 1/2+1/2*(-1-2*a)^(1/2)] [ 1/2-1/2*(-1-2*a)^(1/2)] v = [ -1/2+1/2*(-1-2*a)^(1/2)] [ -1/2-1/2*(-1-2*a)^(1/2)] 对于有些无法求出解析解的非线性方程组,MATLAB只给出一个数值解。这一点可以 从表示解的数字不被方括号括住而确定。例如:键入: [x,y]=solve('sin(x+y)-exp(x)*y=0','x^2-y=2') 得: x = -6.0173272500593065641097297117905 y = 34.208227234306296508646214438330 由于这两个数字没有被[ ]括住,所以它们是数值解。 另外,可利用solve来解线性方程组的通解。例如:键入 P1='2*x1+7*x2+3*x3+x4=6'; P2='3*x1+5*x2+2*x3+2*x4=4'; P3='9*x1+4*x2+x3+7*x4=2'; u=solve(P1,P2,P3,'x1','x2','x3','x4') Warning: 3 equations in 4 variables. u = x1: [1x1 sym] x2: [1x1 sym] x3: [1x1 sym] x4: [1x1 sym] 可以看到:屏幕提示“有3个方程4个变量”,意为解不唯一。(有时会提示解不唯一)且输出的是解的结构形式。为进一步得到解,可输入: u.x1,u.x2,u.x3,u.x4, 得: ans = x1 ans = -5*x1-4*x4 ans = 11*x1+9*x4+2 ans = x4 这样就得到了原方程组的通解。 ⑷ 解符号微分方程 dsolve('eq1','eq2',…) 其中eq表示相互独立的常微分方程、初始条件或 指定的自变量。默认的自变量为t。如果输入的初 始条件少于方程的个数,则在输出结果中出现常数 c1,c2,等字符。关于微分方程的表达式有如下的约 定:字母y表式函数,Dy表示y对t的一阶导数; Dny表示y对t的n阶导数。 例如: 求   的解。可键入:[x,y]=dsolve('Dx=y','Dy=-x') 得 x = cos(t)*C1+sin(t)*C2 y = -sin(t)*C1+cos(t)*C2 dsolve中的输入宗量最多只能有12个,但这并不妨碍解具有多个方程的方程组,因为 可以把多个方程或初始条件定义为一个符号变量进行输入。 例如求  , , f(0)=0 , g(0)=1 的解。可输入指令: P='Df=3*f+4*g,Dg=-4*f+3*g'; v='f(0)=0,g(0)=1'; [f,g]=dsolve(P,v) f = exp(3*t)*sin(4*t) g = exp(3*t)*cos(4*t) 注意:微分方程表达式中字母D必须大写。 求解微分方程   可输入 y=dsolve('D3y=-y','y(0)=1,Dy(0)=0,D2y(0)=0','x') 得: y = (1/3+2/3*exp(1/2*x)*cos(1/2*3^(1/2)*x)*exp(x))/exp(x) 最后看一个解非线性微分方程的例子: dsolve('(Dy)^2+y^2=1','y(0)=0','x') ans = [ sin(x)] [ -sin(x)] 对于无法求出解析解的非线性微分方程,屏幕将提示出错信息。 与数模有关的例 1.曲线拟合 美国人口预测 下表是美国人口统计数据,根据这份资料预测2000年美国人口总数。 年 1790 1800 1810 1820 1830 1840 1850  人口(百万) 3.9 5.3 7.2 9.6 12.9 17.1 23.2  年 1860 1870 1880 1890 1900 1910 1920  人口(百万) 31.4 38.6 50.2 62.9 76.0 92.0 106.5  年 1930 1940 1950 1960 1970 1980 1990  人口(百万) 123.2 131.7 150.7 179.3 204.0 226.5 251.4  Step 1 A=[1790,1800,1810,…; 3.9, 5.3, 7.2,…]’; Step 2 P=polyfit(A(:,1),A(:,2),3) Step 3 px=poly2str(P,'x') Step 4 polyval(P,2000) 如果想了解 fx与数据对x-y的拟和程度,绘出二者的图形最为直观,为此可键入: ft=polyval(P,A(:,1));plot(A(:,1),A(:,2),'bo',A(:,1),ft,'r-')得图形。图中蓝色小圆圈是数据对的图形;而红线是拟合多项式的图形。最后,可与demo_sensus比较。 插值 “线性插值” linear “三次样条插值” spline “三次多项式插值” cubic 对于以上问题,也可以用这三个命令来做。 交通流量问题 下图给出了某城市部分单行街道的交通流量(每小时过车数) x3 100 x6 300 x4 400 200 x2 x7 300 x1 600 x8 500 200 400 300 500 x9 x10 600 700 所给问题满足下列方程组 x1-x3+x4=300 x4+x5=500 x7-x6=200 x1+x2=800 x1+x5=800 x7+x8=1000 x8+x3+x6=1000 (x9=400, x10=600) Step1 A=[0,1,-1,1,0,0,0,0; 0,0,0,1,1,0,0,0; 0,0,0,0,0,-1,1,0; 1,1,0,0,0,0,0,0; 1,0,0,0,1,0,0,0; 0,0,0,0,0,0,1,1; 0,0,1,0,0,1,0,1]; Step 2 b=[300,500,200,800,800,1000,1000]’; Step 3 B1=rank(A);B2=rank([A,b]); Step 4 X=linsolve(A,b) 得特解。 4.线性规划有约束极小问题 用命令x=lp(C,A,b,vlb,vub)。 例:Find x that minimizes f(x)=-5x1-4x2-6x3 subject to x1-x2+x3≦20 3x1+2x2+4x3≦42 3x1+2x2≦30 0≦x1, 0≦x2,0≦x3 First, enter the coefficients: f = [-5; -4; -6] A = [1 -1 1 3 2 4 3 2 0]; b = [20; 42; 30]; lb = zeros(3,1); Next, call a linear programming routine: x= lp(f,A,b,lb); Entering x x = 0.0000 15.0000 3.0000 实际此命令改为: x = linprog(f,A,b,Aeq,beq) x = linprog(f,A,b,Aeq,beq,lb,ub) 对以上的问题可做如下的操作: First, enter the coefficients: f = [-5; -4; -6] A = [1 -1 1 3 2 4 3 2 0]; b = [20; 42; 30]; lb = zeros(3,1); Next, call a linear programming routine: [x,fval,exitflag,output,lambda] = linprog(f,A,b,[],[],lb); Entering x, fval,lambda.ineqlin, and lambda.lower gets x = 0.0000 15.0000 3.0000 fval = -78.0000 和其它信息。 5.非线性规划有约束极小问题 用命令x=constr('f ',x0)。 Examples Find values of x that minimize f(x)=-x1x2x3, starting at the point x = [10; 10; 10] and subject to the constraints 0≤x1+2x2+2x3≤72. -x1-2x2-2x3≤0,x1+2x2+2x3≤72, 第一步:编写M文件 function [f,g]=myfun(x) f=-x(1)*x(2)*x(3); g(1)=-x(1)-2*x(2)-2*x(3); g(2)=x(1)+2*x(2)+2*x(3)-72; 第二步:求解 在MATLAB工作窗中键入 x0=[10,10,10]; x=constr('myfun',x0)即可