Matlab原理与工程应用第二版第五章(函数.ppt

5.1 创建函数,例：x=cos(at)+b y=|x|+c 创建一个comxy函数计算，文件名为comxy.m functionx,y=comxy(a,t,b,c) x=cos(a*t)+b y=abs(x)+c 大写字母开头,5.1.1 函数文件,第五章 函数,1,2,functionx,y=comxy(t,a,b,c) %变量名可不一致 %x=sin(at)+b x=cos(a*t)+b; y=abs(x)+c;,n=3; b=linspace(1,1.4,n); for k=1:3 u,v=comxy(0:pi/4:pi,sqrt(1.8/(1+k)3),b(k),1/0.85) end,调用程序,3,4,functionx,y=comxy(t) %x=sin(at)+b global A B C x=cos(A*t)+B; y=abs(x)+C;,global A B C n=3; C=1/.85; b=linspace(1,1.4,n); for k=1:n A=sqrt(1.8/(1+k)3); B=b(k); m,n=comxy(0:pi/4:pi) end,全局变量的用法，大写，空格（不用逗号）,5,functionx,y=comxy(t,a,b,c) %x=sin(at)+b x=cos(a*t)+b; y=abs(x)+c;,clear all n=3; b=linspace(1,1.4,n); t=0:pi/4:pi; p=zeros(n,length(t); q=zeros(n,length(t); for k=1:3 p(k,:),q(k,:)=comxy(t,sqrt(1.8/(1+k)3),b(k),1/.85) end,6,5.2.2 inline,fox=inline(x.2.*cos(a*x)-b,x,a,b) b=fox(pi/3 pi/3.5,4,1),5.2.3 子函数,7,8,5.3 用户自定义函数、函数句柄和feval,函数句柄是一种引用函数的方法。函数句柄由和函数名构成。 feval就是把已知的数据或符号带入到一个定义好的函数句柄中。,clc clear all syms t %定义符号变量 f=(x,y) x3+y2 k1=feval(f,2,t) % k2=f(2,t) k3=feval(f,1,2) % k4=f(1,2),9,5.4 以数组变量为输入参数的MATLAB函数,5.4.1 多项式拟合-polyfit/polyval,函数polyfit()采用最小二乘法对给定数据进行多项式拟合，其具体使用方法如下： c = polyfit(x,y,n)，采用n次多项式p来拟合数据x和y。 y(x)=c1xn+ c2xn-1 + +cnx + cn+1 (5.1) 求得c后，可以用下列语句求出式5.1 y=polyval(c,xnew),10,11,function nd=data nd=50,0.13;70,0.092;90,0.072;110,0.057;130,0.046;150,0.037;170,0.028;190,0.020;210,0.015;230,0.010;250,0.007;,ncs=data; c=polyfit(ncs(:,1),ncs(:,2),4); r=input(enter r(0r0.2); su=input(enter su(50su250); q=1/(1+polyval(c,su)/sqrt(r); disp(q= num2str(q),函数名以字母开头,p=2 2 3; a=1,0;0,2; m=polyval(p,a),例：求m,12,5.4.2 插值-interp1,一维插值就是对函数y=f(x)进行插值，一维插值的原理如下图所示。,13,函数interp1()实现一维插值，其具体使用方法如下： yi=interp1(x,y,xi)，x，y是已知数据集且具有相同长度的向量； yi = interp1(y,xi)，默认x为1:n，其中n为向量y的长度； yi = interp1(x,y,xi,method)，method用于指定插值的方法。,14,subplot(1,3,1); plot(x,y,ro,xi,yi_nearest,b-); title(最邻近插值); subplot(1,3,2); plot(x,y,ro,xi,yi_linear,b-); title(线性插值); subplot(1,3,3); plot(x,y,ro,xi,yi_spline,b-); title(三次样条插值);,x=0:1.2:10; y=sin(x); xi=0:0.1:10; yi_nearest=interp1(x,y,xi,nearest); yi_linear=interp1(x,y,xi); yi_spline=interp1(x,y,xi,spline); figure; hold on;,15,非线性方程数值求解 函数fzero()求一元函数的零点，其具体使用方法如下： x = fzero(fun,x0,options,p1,p2,)，在x0点附近寻找函数的零点； x = fzero(fun,x0,x1)，在x0,x1区间内寻找函数的零点； x = fzero(fun,x0,options)，用options指定寻找零点的优化器参数； x,fval = fzero(.)，附加自变量为x时的函数值。,16,例：求f(x)=x-10x+2=0在x0=0.5附近的根。 步骤如下： (1) 建立函数文件funx.m。 function fx=funx(x) fx=x-10.x+2; (2) 调用fzero函数求根。 z=fzero(funx,0.5) z = 0.3758 (2) zero=inline(x-10.x+2,x) z=fzero(zero,0.5) (3) z=fzero(inline(x-10.x+2,x) ,0.5),17,对于多值函数，应使用,x1,x2形式，并且准确划定范围，然而，如果f（x1）, f（x2）符号相同，则会得到一个错误。 例：w=fzero(cos,0.6*pi,2*pi)/pi w=fzero(cos,0,2*pi)/pi,18,非线性方程组的求解 对于非线性方程组F(X)=0，用fsolve函数求其数值解。fsolve函数的调用格式为： X=fsolve(funnctionname,X0,option，p1,p2,) 其中X为返回的解，funnctionname是用于定义需求解的非线性方程组的函数文件名，X0是求根过程的初值，option为最优化工具箱的选项设定。最优化工具箱提供了20多个选项，用户可以使用optimset命令将它们显示出来。如果想改变其中某个选项，则可以调用optimset()函数来完成。例如，Display选项决定函数调用时中间结果的显示方式，其中off为不显示，iter表示每步都显示，final只显示最终结果。optimset(Display,off)将设定Display选项为off。,19,求下列非线性方程组在(0.5,0.5) 附近的数值解。 (1) 建立函数文件myfun.m。 function q=myfun(p) x=p(1); y=p(2); q(1)=x-0.6*sin(x)-0.3*cos(y); q(2)=y-0.6*cos(x)+0.3*sin(y); (2) 在给定的初值x0=0.5,y0=0.5下，调用fsolve函数求方程的根。 x=fsolve(myfun,0.5,0.5,optimset(Display,off) x = 0.6354 0.3734,20,将求得的解代回原方程，可以检验结果是否正确，命令如下： q=myfun(x) q = 1.0e-009 * 0.2375 0.2957 可见得到了较高精度的结果。,21,（1）多项式的根 函数roots()求多项式的根，其具体使用方法如下： r = roots(c)，返回多项式c的所有根r。 （2）由根创建多项式 函数poly()实现由根创建多项式，其具体使用方法如下： p = poly(r)，输入r是多项式所有根，返回值为多项式的系数向量； p = poly(A)，输入A是方阵，返回值为A的特征多项式的系数向量。,22,在MATLAB 7语言里，多项式由一个行向量表示，设为p，它的系数按降序排列，使用roots函数可以求出该多项式的根。其使用格式为roots(p)。 p=1 0 3 12 -7 %f(x)=x4+3x2+12x-7 p = 1 0 3 12 -7 roots(p) ans = 0.7876 + 2.4351i 0.7876 - 2.4351i -2.0872 0.5121 ,23,如果已知某个多项式的根，那么，使用poly函数，可以很轻松地产生其对应的多项式。 r=roots(1,-10,35,-50,24) c=poly(r) 使用conv函数对多项式进行乘法运算。 格式为c=conv(a,b)，其中a和b为两个多项式的系数向量，c为相乘所生成的多项式的系数向量。 a=1 2 3 4; b=5 6 7 8; c=conv(a,b),24,8.1.2 数值积分的实现方法 1变步长辛普生法 基于变步长辛普生法，MATLAB给出了quad函数来求定积分。该函数的调用格式为： I,n=quad(fname,a,b,tol,trace) 其中fname是被积函数名。a和b分别是定积分的下限和上限。tol用来控制积分精度，缺省时取tol=0.001。trace控制是否展现积分过程，若取非0则展现积分过程，取0则不展现，缺省时取trace=0。返回参数I即定积分值，n为被积函数的调用次数。,25,求定积分。 (1) 建立被积函数文件fesin.m。 function f=fesin(x) f=exp(-0.5*x).*sin(x+pi/6); (2) 调用数值积分函数quad求定积分。 S,n=quad(fesin,0,3*pi) %n可省略 S = 0.9008 n = 77,26,二重定积分的数值求解 使用MATLAB提供的dblquad函数可以求二重定积分的数值解。该函数的调用格式为： I=dblquad(f,a,b,c,d,tol,trace) 该函数求f(x,y)在a,bc,d区域上的二重定积分。参数tol，trace的用法与函数quad完全相同。,27,计算二重定积分 (1) 建立一个函数文件fxy.m： function f=fxy(x,y) global Ki; Ki=Ki+1; %ki用于统计被积函数的调用次数 f=exp(-x.2/2).*sin(x.2+y); (2) 调用dblquad函数求解。 global Ki; Ki=0; I=dblquad(fxy,-2,2,-1,1) Ki,28,常微分方程（Ordinary differential equations, ODE),初值问题-给出初始值 边值问题-给出边界条件,与初值常微分方程解算有关的指令 ode23 ode45 ode113 ode23t ode15s ode23s ode23tb,29,常用调用格式：t,y=ode45(odefun,tspan,y0) 参数说明： odefun 表示f(t,y)的函数句柄或inline函数，t是标量，y是标量或向量； tspan 如果是二维向量t0,tf，表示自变量初值t0和tf；如果是高维向量t0,t1,tn，则表示输出结点列向量； y0 表示初值向量y0； t 表示结点列向量(t0,t1,tn)T; y 表示数值解矩阵，每一列对应y的一个分量； 完整调用格式： t,y=ode45(odefun,tspan,y0，options,p1,p2,) options 为计算参数设置（如精度要求等），默认表示； p1，p2， 附加传递参数，odefun表示为f(t,y,p1,p2,),30,一.解ODE的基本机理:,2. 把高阶方程转换成一阶微分方程组,1. 列出微分方程,初始条件,令,(5.2),(5.3),(5.4),31,例：著名的Van der Pol方程,令,降为一阶,初始条件,32,3. 根据式(5.3)编写计算导数的M函数文件-ODE文件,把t,Y作为输入宗量，把 作为输出宗量,%M function file name: dYdt.m function Yd = f (t, Y) Yd = f (t,Y) 的展开式,例Van der Pol方程,%M function file name: dYdt.m function Yd=f(t,Y) Yd=zeros(size(Y); Yd(1)=Y(2); Yd(2)=-(Y(1).2-1)*Y(2)-Y(1);,33,4. 使编写好的ODE函数文件和初值 供微分方程解算指令（solver）调用,Solver解算指令的使用格式,输出宗量形式,说明： t0：初始时刻；tN：终点时刻Y0：初值； tol：计算精度,34,例题1：著名的Van der Pol方程,% 主程序 （程序名：VanderPol _ex1.m） t0=0; tN=20; tol=1e-6; Y0=0.25; 0.0; t,Y=ode45(f,t0,tN, Y0,tol); subplot (121), plot (t,Y) subplot (122), plot (Y(:,1), Y(:,2),解法：采用ODE命令,35,Van der Pol方程,% 子程序 （程序名： dYdt.m ） function Ydot = dYdt (t, Y),Ydot=Y(2);-Y(2)*(Y(1)2-1)-Y(1);,或写为,function Ydot = dYdt (t, Y) Ydot=zeros(size(Y); Ydot(1)=Y(2); Ydot(2)=-Y(2)*(Y(1).2-1)-Y(1);,36,37,常微分方程一阶方程组边值问题MATLAB标准形式： 调用格式 Sinit=bvpinit(tinit,yinit) 由在粗略结点tinit的预估解yinit生成粗略解网格sinit Sol=bvp4c(odefun,bcfun,sinit) odefun是微分方程组函数，bcfun为边值条件函数，sol.x为求解结点，sol.y是y(t)的数值解 Sx=deval(sol,ti) 计算由bvp4c得到的解在ti的值,38,例题：求解边值问题 解：首先改写成方程组 边界条件为,39,求解用M函数eg6_5fun.m % M函数eg6_5fun.m clear;close; sinit=bvpinit(0:4,1;0) odefun=inline(y(2);-abs(y(1),t,y); bcfun=inline(ya(1);yb(1)+2,ya,yb); sol=bvp4c(odefun,bcfun,sinit) t=linspace(0,4,101); y=deval(sol,t); plot(t,y(1,:),sol.x,sol.y(1,:),o,sinit.x, sinit.y(1,:),s) legend(解曲线,解点,粗略解),40,最优化是求最优解，也就是在某个区间内有条件约束或者无条件约束地找到函数的最大值或者最小值。MATLAB使用数字方法求函数的最小值。使用迭代算法，也就是有些步骤要重复许多次。现在，假设要求函数f在某个区间内的最小值xmin。 迭代方法需要一个初始估计值x0。从x0开始找到一个更接近xmin的新值x1，这个值的好坏取决于使用的数学方法。直到找到有足够精度的近似值xi才停止迭代，也就是绝对值|xminxi|足够小。 这里提到了标准MATLAB系统的两个最优化命令，fmin命令可以求单变量函数的最小值；fmins命令可以求多变量函数的最小值，同时它还要求有一个初始向量。 在新版本中，fmin和fmins分别被fminbnd和fminsearch取代。 Matlab中没有求函数f的最大值的命令，但可以通过求其相反函数h=f 的最小值间接求得。,41,函数格式： x=fminbnd( fcn , x1 , x2 , options) 求函数在区间(x1，x2 )内取最小值时的x值，采用黄金分割法和二次插值。fcn是目标函数名。如果没有局部最小值，则返回区间内的最小x值。向量options为控制参数，如options(1) = 1，显示中间结果；options(2)表示得到的结果x的精度，缺省为104，options可省略。 x=fmins (fcn, x0, options) 求函数fcn的最小值。由用户自己给出一个初始估计向量x0，优化参数options可省略。,42,例：在区间 0，2 内求函数f(x)=x3-2x-5的最小值。 f=inline(x.3-2*x-5); x = fminbnd(f, 0, 2) x = 0.8165 f(x) ans = -6.0887,符号数学工具箱是操作和解决符号表达式的符号数学工具箱(函数)集合，有复合、简化、微分、积分以及求解代数方程和微分方程的工具。另外还有一些用于线性代数的工具，求解逆、行列式、正则型式的精确结果，找出符号矩阵的特征值而无由数值计算引入的误差。 符号数学工具箱中的工具是建立在功能强大的称作Maple软件的基础上。它最初是由加拿大的滑铁卢（Waterloo）大学开发的。当要求MATLAB进行符号运算时，它就请求Maple去计算并将结果返回到MATLAB命令窗口。因此，在MATLAB中的符号运算是MATLAB处理数字的自然扩展。,44,定义符号变量 参与符号运算的对象可以是符号变量、符号表达式或符号矩阵。符号变量要先定义，后引用。 1. sym函数：每次定义一个符号变量。 x=sym(x) %创建变量x。 x=sym(x, real) %创建实变量x。x=sym(x, unreal)去掉x的实数属性。 a=sym(alpha) %创建变量a，打印为alpha。 2. syms函数 sysms arg1 arg2 argN。一次可定义多个符号变量，使用起来更方便。注意各变量间用空格隔开。,45,符号表达式 符号表达式的书写格式与数值表达式一样，由符号变量、函数、运算符组成，如： syms a b c x f=a*x2 + b*x + c 也可以用f = sym(a*x2 + b*x + c) 把符号表达式ax2+bx+c赋给f，但此时没有定义变量a b c x。 有两种方法创建符号函数，一是用符号表达式，如上面f=a*x2 + b*x + c 定义了符号函数f，从而可以用diff等符号函数对f进行运算。如diff(f)将返回结果2*a*x+b。,46,极限 函数 limit 格式 limit(F,x,a) %计算符号表达式F=F(x)当xa时的极限值。 limit(F,a) %用