1、1234567vpa(exp(sym(1),300)8A=A;1 2 3 1;2;3;4 9【例2-3】试输入复数矩阵需要避免的语句10U S V=svd(X);11【例2-4】用不同的步距生成(0,p)间向量v1=0:0.2:piv2=0:-0.1:piv3=0:piv4=pi:-1:0v5=0:0.2:pi,piv6=pi:-1:0 01213141516矩阵的提取矩阵的提取1、矩阵下标、矩阵下标MATLAB通过确认矩阵下标，可以对矩阵进行插入子块，提取子块和重排子块的操作。1718192021a=1 2 3;3 4 5;m,n=size(a)m=2n=3length(a)ans=3max
2、(size(a)ans=3rank(a)ans=2222324完全下标元素赋值建立三维数组完全下标元素赋值建立三维数组b=1 1;2 2;b(:,:,2)=5；size(b)=2 2 2;bb(:,:,1)=1 1 2 2b(:,:,2)=5 5 5 525ones,zeros,rand和和randn直接创建多维数组的方直接创建多维数组的方法法ones(2,3,3)cat(维，维，p1，p2,)按指定行列数放置模块数组生成多维数组。参数维是指沿着第几维连接模块数组p1和p2等。rand(2,3,3)26cat(维，维，p1，p2,)按指定行列数放置模块数组生成多维数组。参数维是指沿着第几维连接
3、模块数组p1和p2等。a=1 2;3 4;b=5 6;7 8;cat(3,a,b)ans(:,:,1)=1 2 3 4ans(:,:,2)=5 6 7 8cat(1,a,b)ans=1 2 3 4 5 6 7 8cat(2,a,b)ans=1 2 5 6 3 4 7 827 数据结构体将某一类相关的信息纳入一个统一的变量名下数据结构体将某一类相关的信息纳入一个统一的变量名下进行管理。比如要建立学生档案就可以采用数据结构体。结进行管理。比如要建立学生档案就可以采用数据结构体。结构体可以理解成在构体可以理解成在MATLAB语言下的一个小型语言下的一个小型“数据库数据库”例例.考虑建立学生档案结构体
4、，这里应包含下列信息：考虑建立学生档案结构体，这里应包含下列信息：l 编号（用编号（用number表示，而在结构体下表示，而在结构体下number又称为其又称为其成成员变量员变量或域）：数值型。或域）：数值型。l 姓名（姓名（name）：字符串型。）：字符串型。l 身高（身高（height）:数值型数值型l 考试成绩（考试成绩（test）：矩阵，其第）：矩阵，其第i行为第行为第i次考试的成绩，第次考试的成绩，第j列为第列为第j门考试的成绩门考试的成绩结构体的建立、修改及应用结构体的建立、修改及应用28建立结构体建立结构体student_rec student_rec.number=1;stud
5、ent_rec.name=张三张三;student_rec.height=180;student_rec.test=100,80,75;77,60,92;67,28,90;100,89,78;显示结构体的内容显示结构体的内容 student_recstudent_rec=number:1 name:张三张三 height:180 test:4x3 double29和和C语言不同，引用结构体中成员变量时直接引用即可，语言不同，引用结构体中成员变量时直接引用即可，而无需用而无需用-算符来引用。算符来引用。获得该结构体中的获得该结构体中的test成员变量成员变量 student_rec.testan
6、s=100 80 75 77 60 92 67 28 90 100 89 7830建立一组建立一组100个学生（假设个学生（假设2个年级，每个年级个年级，每个年级50人）的新的人）的新的结构体变量结构体变量b b(50,2)=struct(student_rec)%构造构造502结构体矩阵结构体矩阵b=50 x2 struct array with fields:number name height test二年级第二年级第43个学生的有关信息的填写个学生的有关信息的填写 b(43,2).number=50+43;b(43,2).name=李四李四;b(43,2).height=178;b(4
7、3,2).test=83 80 78;97 80 72;69 88 80;87 99 100;31结构体的添加，修改结构体的添加，修改给给b变量加一个成员变量变量加一个成员变量体重（体重（weight）,则可以将结构则可以将结构体中任意一个变量后面加一个体中任意一个变量后面加一个weight成员变量，即成员变量，即 b(1,2).weight=90b=50 x2 struct array with fields:number name hight test weight32成员变量的消除成员变量的消除 rmfield()消除消除b结构体中结构体中weight成员变量成员变量 b=rmfield
8、(b,weight)b=50 x2 struct array with fields:number name hight test3334用单元数据结构来构造某个学生档案用单元数据结构来构造某个学生档案 B=1,张三张三,180,100 80 75;77 60 92;67 28 90;100 89 78B=1 张三张三 180 4x3 double size(B)ans=1 4 B4ans=100 80 75 77 60 92 67 28 90 100 89 78从显示中不能直从显示中不能直接看出此单元数接看出此单元数据的第据的第4个元素代个元素代表的矩阵具体内表的矩阵具体内容是什么容是什么不
9、能用圆括号去不能用圆括号去访问单元数组访问单元数组35使用使用 celldisp()函数显示整个单元函数显示整个单元 celldisp(B)%显示整个单元变量显示整个单元变量 B1=1 B2=张三张三 B3=180 B4=100 80 75 77 60 92 67 28 90 100 89 7836消除单元元素消除单元元素 删除删除B单元变量的第三个单元元素单元变量的第三个单元元素 B(3)=B=1 张三张三 4x3 double这里用的是这里用的是B（3）而）而不是不是B3，后者只能，后者只能将第将第3单元置成空矩阵，单元置成空矩阵，而不能消除它而不能消除它3738复数不变复数求共轭3940
10、4142点运算在点运算在MATLAB中起着非常重要的作用，例如中起着非常重要的作用，例如当当x是一个向量时，则求取数值是一个向量时，则求取数值 ，不能直接写，不能直接写成成 而必须写成而必须写成 。在进行矩阵的点运算。在进行矩阵的点运算时，同样要求运算的两个矩阵的维数一致，或其时，同样要求运算的两个矩阵的维数一致，或其中一个变量为标量。其实一些特殊的函数，如中一个变量为标量。其实一些特殊的函数，如sin()也是由点运算的形式进行的，因为它要对矩阵的也是由点运算的形式进行的，因为它要对矩阵的每个元素求正弦值。每个元素求正弦值。5x5x5.x43+加加.数组幂数组幂-减减左除或反斜杠左除或反斜杠*
11、矩阵乘法矩阵乘法/右除或斜杠右除或斜杠.*数组乘法数组乘法./数组除数组除矩阵幂矩阵幂：冒号冒号4445Matlab还提供一些特殊的函数，在编程中很实用还提供一些特殊的函数，在编程中很实用n find()函数函数 可以查询出满足某种关系数组的下标可以查询出满足某种关系数组的下标 C=1 1 1 1;0 1 1 1;find(C=1)ans=1 3 4 5 6 7 8 i,j=find(C=1);i,jans=1 1 1 2 2 2 1 3 2 3 1 4 2 446n查询函数查询函数all()和和any()C=1 1 1 1;0 1 1 1;all(C=1)ans=0 1 1 1 C=1 1
12、1 1;0 1 1 1;any(C=1)ans=1 1 1 1all()和和any()函数可以互相嵌套函数可以互相嵌套例如，判定一个矩阵例如，判定一个矩阵C是否元素均非零可以使用是否元素均非零可以使用all(all(C)all(C(:)47矩阵元素的数据变换矩阵元素的数据变换对由小数构成的矩阵对由小数构成的矩阵A来说，如果想对它取整，有来说，如果想对它取整，有n floor(A)将将A中的元素按负无穷方向取整中的元素按负无穷方向取整nceil(A)将将A中元素按正无穷方向取整中元素按正无穷方向取整nround(A)将将A中元素按最近的整数取整，即中元素按最近的整数取整，即 四舍五入四舍五入nf
13、ix(A)将将A中元素按离中元素按离0近的方向取整近的方向取整n求取矩阵元素的余数求取矩阵元素的余数 C=rem(A,x)n将小数表示成两个整数的除式形式将小数表示成两个整数的除式形式n,d=rat(A),A=n./d48 floor(A)%向向-Inf方向取整方向取整ans=-3 2 1 1 -1 0 -1 -1 -2 A=-2.5+5*rand(3)%生成生成-2.5 2.5区间上均匀分布的随机数区间上均匀分布的随机数A=-2.42363036485482 2.15907289230832 1.73110708912162 1.23392838282215 -0.1700282916228
14、8 0.12576248152586 -0.27451783856027 -0.40675266136247 -1.48676321174806 ceil(A)%向向+Inf方向取整方向取整ans=-2 3 2 2 0 1 0 0 -1 round(A)%取最近的整取最近的整数数ans=-2 2 2 1 0 0 0 0 -1 fix(A)%向向0方向取整方向取整ans=-2 2 1 1 0 0 0 0 -149(2)5051p 最小二乘曲线拟合技术最小二乘曲线拟合技术p 代数方程的求解代数方程的求解p 微分方程的求解微分方程的求解p 无约束最优化问题求解无约束最优化问题求解p 二次型规划问题求
15、解二次型规划问题求解p 一般非线性规划问题求解一般非线性规划问题求解525354555657多项式拟合函数585960二维数据的插值拟合二维数据的插值拟合二维插值函数二维插值函数61x,y=meshgrid(-3:.6:3,-2:.4:2);z=(x.2-2*x).*exp(-x.2-y.2-x.*y);mesh(x,y,z),axis(-3,3,-2,2,-0.7,1.5)二元函数二元函数222(,)(2)xyxyzf x yxx e绘制已知数据的网格图Exp3_5062默认插值算法默认插值算法立方和样条插值算法立方和样条插值算法x1,y1=meshgrid(-3:.2:3,-2:.2:2)
16、;z1=interp2(x,y,z,x1,y1);mesh(x1,y1,z1),axis(-3,3,-2,2,-0.7,1.5)z1=interp2(x,y,z,x1,y1,cubic);z2=interp2(x,y,z,x1,y1,spline);mesh(x1,y1,z1),axis(-3,3,-2,2,-0.7,1.5)figure;mesh(x1,y1,z2),axis(-3,3,-2,2,-0.7,1.5)63MATLAB提供的二维插值函数能比较好地进行提供的二维插值函数能比较好地进行二维插值运算，它只能处理以网格形式给出的二维插值运算，它只能处理以网格形式给出的数据，如果已知数据不
17、是以网格形式给出的，数据，如果已知数据不是以网格形式给出的，则用该函数是无能为力的。在实际应用中大部则用该函数是无能为力的。在实际应用中大部分问题都是以实测的分问题都是以实测的 给出的，所给出的，所以不能直接使用该函数进行二维插值。以不能直接使用该函数进行二维插值。(,)iiix y z6465根据函数原型建立根据函数原型建立MATLAB函数函数求待定系数向量求待定系数向量待定系数待定系数向量向量a目标函数目标函数值值原型函数原型函数名名优化初值优化初值原始输入原始输入输出数据输出数据向量向量Exp3_5166【例2-38】用最小二乘拟合方法拟合函数的待定系数67examp2_38指定函数最指
18、定函数最小小/最大值最大值686970322121xx138 =a x =b例:x1+2x2=8 2x1+3x2=13717221xx321 =433221 a x =b7374432321321xxx21=a x =b7500)(),(yxyyxfdtdy76111(,)nnnnnnnnyyyydyf xydtxxh7778采取措施采取措施n选择适当的步长选择适当的步长 采用像采用像Euler法这样简单的算法时，法这样简单的算法时，应适当地选择步长，既不能太大，又不能太小；应适当地选择步长，既不能太大，又不能太小；n改进近似算法精度改进近似算法精度 由于由于Euler算法只是将原积分问算法只
19、是将原积分问题进行梯形的近似，其近似精度很低，不能有效逼近题进行梯形的近似，其近似精度很低，不能有效逼近原始问题。可以采用各种更精确的插值方法来取代原始问题。可以采用各种更精确的插值方法来取代Euler公式，从而改进运算精度。如公式，从而改进运算精度。如Runge-Kutta法、法、Adams法等。法等。n采用变步长方法采用变步长方法 前面提到前面提到“适当适当”地选择步长本身地选择步长本身就是模糊的概念，步长的选取决定于经验。事实上，就是模糊的概念，步长的选取决定于经验。事实上，很多种方法都允许变补偿的求解，如果误差较小时，很多种方法都允许变补偿的求解，如果误差较小时，可自动地增大步长，而误
20、差较大时，再自动的减小步可自动地增大步长，而误差较大时，再自动的减小步长，从而精确、有效地求解给出的常微分方程初值问长，从而精确、有效地求解给出的常微分方程初值问题。题。792112121kkyynnk1=hf(xn,yn)k2=hf(xn+h,yn+k1)8081四阶五级四阶五级Runge-Kutta变步长算法变步长算法德国学者德国学者Felhberg对传统的对传统的Runge-Kutta法进行了改进法进行了改进误差向量，用它的大小来误差向量，用它的大小来变换计算步长变换计算步长82通用调用格式：83【例2-34】微分方程组84建立函数建立函数在命令窗口输入在命令窗口输入看轨迹动画效果看轨迹
21、动画效果Example_rossler85附加参数的 M-函数优点：若想改变参数，没必要修改优点：若想改变参数，没必要修改M函数，只需函数，只需要在求解该方程时，将参数带入即可，这样会很要在求解该方程时，将参数带入即可，这样会很方便方便8687 ode45(),ode15s()88【例2-35】方程89051015202530-3-2-10123-2.5-2-1.5-1-0.500.511.522.5-3-2-1012390function y=vdp_eq(t,x,flag,mu)y=x(2);-mu*(x(1).2-1).*x(2)-x(1);建立方程建立方程h_opt=odeset;x0
22、=-0.2;-0.7;t_final=20;mu=1;t1,y1=ode45(vdp_eq,0,t_final,x0,h_opt,mu);mu=2;t2,y2=ode45(vdp_eq,0,t_final,x0,h_opt,mu);plot(t1,y1,t2,y2,:)figure;plot(y1(:,1),y1(:,2),y2(:,1),y2(:,2),:)解方程（参数传递）解方程（参数传递）Example3_3091929394function dx=apolloeq(t,x)mu=1/82.45;mu1=1-mu;r1=sqrt(x(1)+mu)2+x(3)2);r2=sqrt(x(1)
23、-u1)2+x(3)2);dx=x(2);2*x(4)+x(1)-mu1*(x(1)+mu)/r13-mu*(x(1)-mu1)/r23;x(4);-2*x(2)+x(3)-mu1*x(3)/r13-mu*x(3)/r23;x0=1.2;0;0;-1.04935751;t,y=ode45(apolloeq,0,20,x0);length(t),plot(y(:,1),y(:,3)example3_3595函数句柄的创建函数句柄的创建函数句柄是包含了函数的路径、函数名、类型等信息。创建使用符号或者str2func实现特点：特点：（1）在更大范围内调用函数。包含了函数文件的路径和函数类型。无）在更
24、大范围内调用函数。包含了函数文件的路径和函数类型。无论函数所在文件是否在搜索路径上，是否当前路径，是否是子函数或论函数所在文件是否在搜索路径上，是否当前路径，是否是子函数或私有函数等，只要函数句柄存在，函数都能执行；私有函数等，只要函数句柄存在，函数都能执行；（2）提高函数调用速度。不使用函数句柄时，对函数的每次调用都要）提高函数调用速度。不使用函数句柄时，对函数的每次调用都要为该函数进行全面的路径搜索，影响速度；为该函数进行全面的路径搜索，影响速度；（3）使函数调用像使用变量一样方便、简单；）使函数调用像使用变量一样方便、简单；（4）可获得同名重载函数的位置、类型、信息）可获得同名重载函数的
25、位置、类型、信息96Exp_3.29Chenshundun.mdlChenchaos.mSynchronization1.m97方程描述方程描述function dx=c3ximp(t,x)A=sin(x(1)cos(x(2);-cos(x(2)sin(x(1);B=1-x(1);-x(2);dx=inv(A)*B;Exp3_3198t,x=ode45(c3ximp,0,10,0;0);plot(t,x)求解语句求解语句99p在用微分方程描述的一个变化过程中，若往往又包含着多个相互作用但变化速度相差十分悬殊的子过程，这样一类过程就认为具有“刚性”。描述这类过程的微分方程初值问题称为“刚性问题”
26、。p例 如，宇航飞行器自动控制系统一般包含两个相互作用但效应速度相差十分悬殊的子系统，一个是控制飞行器质心运动的系统，当飞行器速度较大时，质心运动惯性较大，因而相对来说变化缓慢变化缓慢；另一个是控制飞行器运动姿态的系统，由于惯性小，相对来说变化很快变化很快，因而整 个系统就是一个刚性系统。刚性问题描述刚性问题描述100刚性方程的刚性方程的MATLAB求解求解1000时的时的Van der pol方程是典型的刚性方程方程是典型的刚性方程 examp3_30就不行就不行h_opt=odeset;x0=2;0;t_final=3000;tic,mu=1000;t,y=ode15s(vdp_eq1,0
27、,t_final,x0,h_opt,mu);tocplot(t,y(:,1);figure;plot(t,y(:,2)应该采用应该采用examp3_32(x)101t,y=ode45(c3xstf1,0,1,1;0;-1);x1=exp(-2*t);x2=exp(-40*t).*cos(40*t);x3=exp(-40*t).*sin(40*t);y1=0.5*x1+0.5*x2+0.5*x3,0.5*x1-0.5*x2-0.5*x3,-x2+x3;plot(t,y,r,t,y1,b)function dy=c3xstf1(t,x)dy=-21,19,-20;19,-21,20;40,-40,
28、-40*x;因为在Matlab下采用变步长计算，可以自动按照精度修正，所以感受不到它是一个刚性问题Example3_33102MATLAB的最优化工具箱中提供了多元方程的求解函数的最优化工具箱中提供了多元方程的求解函数fsolve()非线性方程的求解非线性方程的求解返回的解返回的解原函数在原函数在x点处的值点处的值函数返回的条函数返回的条件件1，0解的附加解的附加信息信息求解问题的数求解问题的数学描述学描述自变量起自变量起始搜索点始搜索点优化工具箱的优化工具箱的选项设定选项设定103optimset 优化工具箱常用选项优化工具箱常用选项104function q=my2deq(x)q=zero
29、s(2,1);q(1)=x(1)*x(1)+x(2)*x(2)-1;q(2)=0.75*x(1)3-x(2)+0.9;OPT=optimset;OPT.LargeScale=off;x,Y,c,d=fsolve(my2deq,1;2,OPT),105 Optimization terminated successfully:First-order optimality is less than options.TolFun.x=0.3570 0.9341Y=1.0e-009*0.1215 0.0964c=1d=iterations:7 funcCount:21 algorithm:trust-
30、region dogleg firstorderopt:1.3061e-010初值的选择有时对整个问题的求解有很大的影响，在某初值的选择有时对整个问题的求解有很大的影响，在某些初值下甚至无法搜索到方程的解些初值下甚至无法搜索到方程的解106MATLAB最优化工具箱中的函数只能解决实数解问最优化工具箱中的函数只能解决实数解问题，再配以符号运算工具箱则可以解出所有的解题，再配以符号运算工具箱则可以解出所有的解 syms x y;X=solve(x2+y2-1=0,0.75*x3-y+0.9=0)107解析求解方法并不是万能的，因为这里的例子最终解析求解方法并不是万能的，因为这里的例子最终可以转换为
31、一元高次代数方程所以能用它求解，但可以转换为一元高次代数方程所以能用它求解，但更一般的方程是不能解出的。更一般的方程是不能解出的。X=solve(x2*sin(x)=x+1,y=x*cos(0.1*x2+3*x)+0.5);X.x;X.y 108 ezplot(x2*sin(y)-x-1);hold on ezplot(-y+x*cos(0.1*x2+3*x)+0.5)方程求根过程中不能过分依赖数值方法和解析方方程求根过程中不能过分依赖数值方法和解析方法，有时还应该适当考虑图解的方法法，有时还应该适当考虑图解的方法109110目标函数的目标函数的MATLAB表示表示function y=opt
32、_fun0(x)y=3*x(1)*x(1)-2*x(1)*x(2)+x(2)*x(2)+4*x(1)+3*x(2);求解最优化问题求解最优化问题format long;x,f_opt,c,d=fminsearch(opt_fun0,-1.2,1);x,f_opt,kk=d.funcCountExp_3.45111线性规划问题的数学描述为线性规划问题的数学描述为.minx s t AxBf x调用格式调用格式eqeqA xB线性等式约束线性等式约束线性不等式约束线性不等式约束AxB上界向量与下界向量上界向量与下界向量mMxxx112将之转换为最小化问题将之转换为最小化问题1,2,1,3f线性等式
33、约束线性等式约束1,1,3,16x线性不等式约束线性不等式约束0211301614x 上界向量与下界向量上界向量与下界向量0ix Exp_3.46113f=1,-2,1,-3;Aeq=1,1,3,1;Beq=6;A=0,-2,1,1;0,-1,6,-1;B=3;4;LB=zeros(4,1);UB=;x,f_opt=linprog(f,A,B,Aeq,Beq,LB,UB)Optimization terminated successfully.x=0.00000000000000 1.00000000000000 0.00000000000000 4.99999999999999f_opt=-
34、16.99999999999998114二次型规划问题二次型规划问题二次型规划的数学表示为二次型规划的数学表示为二次型规划问题函数二次型规划问题函数 quadprog()11511115332110 x0ix Exp_3.47116f=-2,-4,-6,-8;H=diag(2,2,2,2);OPT=optimset;OPT.LargeScale=off;%不使用大规模问题求解不使用大规模问题求解A=1,1,1,1;3,3,2,1;B=5;10;Aeq=;Beq=;LB=zeros(4,1);x,f_opt=quadprog(H,f,A,B,Aeq,Beq,LB,OPT)Optimization
35、 terminated successfully.x=0.00000000000000 0.66666666666667 1.66666666666667 2.66666666666667f_opt=-23.66666666666666别忘记补上已经除别忘记补上已经除去了的常数项去了的常数项117118MATLAB提供的求解各种约束下的最优化问题函数提供的求解各种约束下的最优化问题函数 fmincon()目标函数目标函数的的M函数函数初始搜索初始搜索点点不等式约不等式约束束等式约束等式约束上下界上下界非线性约束非线性约束的的M函数函数控制选项控制选项119写出目标函数写出目标函数functio
36、n y=opt_fun1(x)y=1000-x(1)*x(1)-2*x(2)*x(2)-x(3)*x(3)-x(1)*x(2)-x(1)*x(3);写出非线性约束写出非线性约束function c,ceq=opt_con1(x)ceq=x(1)*x(1)+x(2)*x(2)+x(3)*x(3)-25;8*x(1)+14*x(2)+7*x(3)-56;c=;不等式约束不等式约束非线性等非线性等式约束式约束Exp_3.48120OPT=optimset;OPT.LargeScale=off;x0=1;1;1;LB=0;0;0;UB=Inf;Inf;Inf;A=;B=;Ae=;Be=;x,f_opt
37、,c,d=fmincon(opt_fun1,x0,A,B,Ae,Be,LB,UB,opt_con1,OPT);x,f_opt,kk=d.funcCount Optimization terminated successfully:Magnitude of directional derivative in search direction less than 2*options.TolFun and maximum constraint violation is less than options.TolConActive Constraints:1 2x=3.51212879449002 0
38、.21698734902870 3.55216382252544f_opt=9.617151721291774e+002kk=61121考虑到第二个约束条件实际上是线性等式约束考虑到第二个约束条件实际上是线性等式约束function c,ceq=opt_con2(x)ceq=x(1)*x(1)+x(2)*x(2)+x(3)*x(3)-25;c=;非线性约束函数可以写为非线性约束函数可以写为x0=1;1;1;Ae=8,14,7;Be=56;x,f_opt,c,d=fmincon(opt_fun1,x0,A,B,Ae,Be,LB,UB,opt_con2,OPT);x,f_opt,kk=d.funcCountFmincon()函数调用函数调用122123
