MATLAB在数学建模中的应用ppt课件

1、MATLAB软件及其在 数学建模中的应用,1,求解结果 发现规律 模型验证 讨论分析,计算在数学建模中的作用,2,数学建模中的计算,问题的分析,修正模型,粗假设,修正算法,结果分析,讨论推广,修正假设,粗模型,粗算法,发现问题,发现规律,模型验证,3,主要内容,Matlab软件简介 数学建模Matlab算法,4,MATLAB简介,MATLAB是MATrix LABoratory 的缩写，是由美国MathWorks公司开发的工程计算软件，迄今MATLAB已推出了6.5版. 1984年MathWorks公司正式将MATLAB推向市场，从这时起，MATLAB的内核采用C语言编写，而且除原有的数值计算

2、能力外，还新增了数据图视功能.在国际学术界，MATLAB已经被确认为准确、可靠的科学计算标准软件.在设计研究单位和工业部门，MATLAB被认作进行高效研究、开发的首选软件工具.,5,MATLAB的发展 1984年，MATLAB第1版(DOS版) 1992年，MATLAB4.0版 1994年，MATLAB 4.2版 1997年，MATLAB 5.0版 1999年，MATLAB 5.3版 2000年，MATLAB 6.0版 2001年，MATLAB 6.1版 2002年，MATLAB 6.5版 2004年，MATLAB 7.0版,6,MATLAB的功能,MATLAB产品组是从支持概念设计、算法开发

3、、建模仿真， 到实时实现的集成环境，可用来进行： 数据分析 数值与符号计算 工程与科学绘图 控制系统设计 数字图像信号处理 建模、仿真、原型开发 财务工程、应用开发、图形用户界面设计,功能强大,7,MATLAB语言特点,编程效率高，允许用数学的语言来编写程序 用户使用方便，把程序的编辑、编译、连接和执行融为一体 高效方便的矩阵和数组运算 语句简单，内涵丰富 扩充能力强，交互性，开放性 方便的绘图功能 该软件由c语言编写，移植性好,语言简洁,8,学习该软件的必要性：目前，MATLAB软件不仅走入企业、公司和科研机构，而且在高等院校也是从大学生到博士生都必须掌握的一项基本技能，是必不可少的计算工具

4、，。 MATLAB功能：数值计算、符号运算和图形处理。,9,学习它的意义：随着计算机科学和计算软件的发展，数学系学生必须掌握一门好的计算软件。这是我们就业、继续身造或做科研工作所要用到的。是当代大学生必备的一项技能。,10,其它计算软件：MATHEMATIC（数学分析问题的计算）；LINGO（规划问题的计算）。可以说一个人掌握了一门计算软件，再学习其它计算软件就很容易。,11,MATLAB的环境,菜单项； 工具栏； 【Command Window】命令窗口； 【Launch Pad】分类帮助窗口； 【Workspace】工作区窗口； 【Command History】指令历史记录窗口； 【Cu

5、rrent Directory】当前目录选择窗口；,12,MATLAB操作窗口,接受命令的窗口,13,MATLAB在微积分中的应用,1、求函数值,例1 在命令窗口中键入表达式 并求 时的函数值。, x=2,y=4 z=x2+exp(x+y)-y*log(x)-3,x = 2 y = 4 z = 401.6562,命令窗口显示结果：,14,例2 用循环语句编写M文件计算ex的值，其中x,n为输入 变量，ex的近似表达式为,function y=e(x,n) y=1;s=1; for i=1:n s=s*i; y=y+xi/s; end y, y=e(1,100) ans = y y = 2.71

6、83,调用函数 M文件,15,MATLAB在微积分中的应用,2、求极限,例3 求极限, syms n; limit(sqrt(n+sqrt(n)-sqrt(n),n,inf),ans = 1/2,LIMIT Limit of an expression. LIMIT(F,x,a) takes the limit of the symbolic expression F as x - a. LIMIT(F,x,a,right) or LIMIT(F,x,a,left) specify the direction of a one-sided limit.,定义符号变量,16,MATLAB在微积分

7、中的应用,3、求导数, syms x y=10 x+x10+log(x) y = x10+10 x+log(x) diff(y),ans = 10*x9+10 x*log(10)+1/x,定义X为符号变量,求,Difference：差分 Differential：微分的,17, syms x; y=log(1+x); a=diff(y,x,2) a = -1/(1+x)2 x=1;eval(a) ans = -0.2500,求,求,将符号表达式 转换成数值表达式,18,例6 设,，求, syms x y; z=exp(2*x)*(x+y2+2*y); a=diff(z,x) b=diff(z,

8、y) c=diff(z,x,2) d=diff(z,y,2) e=diff(a,y),19,a =2*exp(2*x)*(x+y2+2*y)+exp(2*x) b =exp(2*x)*(2*y+2) c =4*exp(2*x)*(x+y2+2*y)+4*exp(2*x) d =2*exp(2*x) e =2*exp(2*x)*(2*y+2),20,MATLAB在微积分中的应用,4、求极值和零点, fzero(3*x5-x4+2*x3+x2+3,0),ans = -0.8952,起始点,函数,命令函数, fminbnd(3*x5-x4+2*x3+x2+3,-1,2) ans = -1.1791e

9、-005,21,MATLAB在微积分中的应用,4、求极值和零点, X,FVAL= FMINSEARCH(x(1)2+2.5*sin(x(2)- x(3)*x(1)*x(2)2,1 -1 0),X = 0.0010 -1.5708 0.0008 FVAL =-2.5000,22,MATLAB在微积分中的应用,5、求积分,例9 求不定积分, int(cos(2*x)*cos(3*x),ans =1/2*sin(x)+1/10*sin(5*x),例10 求定积分,Integrate：积分, eval(int(x2*log(x),1,exp(1) ans = 4.5746, x=1:0.01:exp(

10、1); y=x.2.*log(x); trapz(x,y) ans = 4.5137,23,例10 求定积分, int(exp(-x2/2),0,1) ans = 1/2*erf(1/2*2(1/2)*2(1/2)*pi(1/2), x=0:0.01:1; y=exp(-x.2/2); trapz(x,y) ans = 0.8556, y=exp(-x.2/2); quadl(y,0,1) ans = 0.8556,变步长 数值积分,梯形法数值积分,24,MATLAB在微积分中的应用,5、求积分,例11 求二重积分, syms x y; f=y2/x2; int(int(f,x,1/2,2),

11、y,1,2) ans =7/2,符号积分, f=(y.2)./(x.2); dblquad(f,1/2,2,1,2) ans = 3.5000,数值计算,25,MATLAB在微积分中的应用,6、解微分方程,例12 计算初值问题：, dsolve(Dy=x+y,y(0)=1,x),ans =-x-1+2*exp(x),一定要大写,26,MATLAB在微积分中的应用,7、级数问题,例13 求函数 的泰勒展开式，并计算该 函数在x=3.42时的近似值。, syms x; taylor(sin(x)/x,x,10),ans = 1-1/6*x2+1/120*x4-1/5040*x6+1/362880*

12、x8, x=3.42; eval(ans) ans = -0.0753,27,MATLAB在线性代数中的应用,1、矩阵的基本运算,例1 已知, a=4 -2 2;-3 0 5;1 5 3; b=1 3 4;-2 0 -3;2 -1 1; a*b,=AB,28,MATLAB在线性代数中的应用,1、矩阵的基本运算,例1 已知, inv(a) ans = 0.1582 -0.1013 0.0633 -0.0886 -0.0633 0.1646 0.0949 0.1392 0.0380,29,MATLAB在线性代数中的应用,1、矩阵的基本运算,例1 已知, rank(a) ans = 3,30,MAT

13、LAB在线性代数中的应用,1、矩阵的基本运算,例1 已知, a/b ans = 0 0 2.0000 -2.7143 -8.0000 -8.1429 2.4286 3.0000 2.2857,31,MATLAB在线性代数中的应用,1、矩阵的基本运算,例1 已知, ab ans = 0.4873 0.4114 1.0000 0.3671 -0.4304 0 -0.1076 0.2468 0,32,2、解线性方程组, a=1 -1 4 -2;1 -1 -1 2;3 1 7 -2;1 -3 -12 6; rref(a),ans =,1 0 0 0 0 1 0 0 0 0 1 0 0 0 0 1,将矩

14、阵A化为最简阶梯形,R（A）=4=n； 所以方程组只有零解。,RREF Reduced row echelon form,33,2、解线性方程组,34,求齐次方程组 的基础解系, a=2 3 1;1 -2 4;3 8 -2;4 -1 9; b=4;-5;13;-6; c=null(a,r) c = -2 1 1,求非齐次方程组 的一个特解, l u=lu(a); x0=u(lb) x0 = -3124/135 3529/270 2989/270,所以方程组的一般解为,35,3、将矩阵对角化, a=-1 2 0;-2 3 0;3 0 2; v,d=eig(a) v = 0 379/1257 37

15、9/1257 0 379/1257 379/1257 1 -379/419 -379/419 d =2 0 0 0 1 0 0 0 1,A的特征值为2，1，1,36,4、用正交变换化二次型为标准形, a=1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1; format u t=schur(a),u =0.0846 0.4928 0.7071 0.5000 0.0846 0.4928 -0.7071 0.5000 -0.7815 -0.3732 0 0.5000 0.6124 -0.6124 0 0.5000 t = -0.0000 0 0 0 0 -0.0000 0 0 0 0

16、0 0 0 0 0 4.0000,37, a=1 1 1 1；1 1 1 1；1 1 1 1；1 1 1 1; format rat u t=schur(a),u = 596/7049 1095/2222 985/1393 1/2 596/7049 1095/2222 -985/1393 1/2 -1198/1533 -789/2114 0 1/2 1079/1762 -1079/1762 0 1/2 t = * 0 0 0 0 * 0 0 “*”表示 0 0 0 0 近似于零 0 0 0 4,FORMAT RAT Approximation by ratio of small integer

17、s.,38,4、用正交变换化二次型为标准形,结论：作正交变换,则有,39,上机实验题 一、基础型实验 1、计算下列极限,40,2、计算下列导数 (1) (2) (3) (4),41,实验练习,一. 输入A=1,1,1;1,2,3;1,3,6,B=8,1,6;3,5,7;4,9,2,u=3;1;4, 1. A+B; 2. A-B; 3. A*B; 4. A*u; 5. 2A-3B; 6. A2+B2; 7. AB-BA。,二. 求下列矩阵的逆阵并求其行列式的值,1. A=1,3,3;1,4,3;1,3,4; 2. A=1,2,3;2,2,1;3,4,3; 3. A=1,1,1,1;1,1,-1,

18、-1;1,-1,1,-1;1,-1,-1,1; 4. A=1,1,0,0;1,2,0,0;3,7,2,3;2,5,1,2。,三. 解矩阵方程,1.A=2,5;1,3,B=4,-6;2,1,AX=B; 2.A=2,1,-1;2,1,0;1,-1,1,B=1,-1,3;4,3,2;1,-2,5,XA=B; 3.A=1,4;-1,2,B=2,0;-1,1,C=3,1;0,-1,AXB=C; 4.A=0,1,0;1,0,0;0,0,1,B=1,0,0;0,0,1;0,1,0, C=1,-4,3;2,0,-1;1,-2,0,AXB=C.,42,四. 将下列矩阵化为阶梯矩阵,1.A=1,-2,0;-1,1

19、,1;1,3,2; 2.A=0,1;1,0;0,-1;,3.A=1,2,3,4;0,1,2,3;0,0,1,2;0,0,0,1;,4.A=2,1,0,0;3,2,0,0;1,1,3,4;2,-1,2,3.,五.求下列矩阵的秩 1.A=-5,6,-3;3,1,11;4,-2,8; 2.A=1,-2,3,-1;3,-1,5,-3;2,1,2,-2; 3.A=3,1,0,2;1,-1,2,-1;1,3,-4,4; 4.A=1,4,-1,2,2;2,-2,1,1,0;-2,-1,3,2,0.,43,附录：MATLAB软件中部分常用函数表,44,作为一个功能强大的工具软件，Matlab具有很强的图形处理

20、功能，提供了大量的二维、三维图形函数。由于系统采用面向对象的技术和丰富的矩阵运算，所以在图形处理方面即方便又高效。,45,一、 plot数据点绘图命令 命令格式：plot(x,y) 其中x和y为坐标向量 命令功能：以向量x、y为轴，绘制曲线。 【例1】 在区间0X2内，绘制正弦曲线Y=sin(x)，其程序为： x=0:pi/100:2*pi; y=sin(x); plot(x,y),1.二维图形,46,【例2】同时绘制正、余弦两条曲线y1=sin(x)和y2=cos(x)，其程序为： x=0:pi/100:2*pi; y1=sin(x); y2=cos(x); plot(x,y1,x,y2)

21、plot函数还可以为plot(x,y1,x,y2，x,y3，)形式，其功能是以公共向量x为X轴，分别以y1，y2，y3，为Y轴，在同一幅图内绘制出多条曲线。,47,（一）线型与颜色 格式：plot(x,y1,cs,.) 其中c表示颜色， s表示线型。,【例3】 用不同线型和颜色重新绘制例2图形，其程序为： x=0:pi/100:2*pi; y1=sin(x); y2=cos(x); plot(x,y1,go,x,y2,b-.) 其中参数go和b-.表示图形的颜色和线型。g表示绿色，o表示图形线型为圆圈；b表示蓝色，-.表示图形线型为点划线。,48,绘图基本线型和颜色,49,（二）图形标记 在绘

22、制图形的同时，可以对图形加上一些说明，如图形名称、图形某一部分的含义、坐标说明等，将这些操作称为添加图形标记。 title(加图形标题); xlabel(加X轴标记); ylabel(加Y轴标记); text(X,Y,添加文本);,50,（三）设定坐标轴 用户若对坐标系统不满意，可利用axis命令对其重新设定。 axis(xmin xmax ymin ymax) 设定最大和最小值 axis （auto） 将坐标系统返回到自动缺省状态 axis （square） 将当前图形设置为方形 axis （equal） 两个坐标因子设成相等 axis （off） 关闭坐标系统 axis （on） 显示坐标

23、系统,51,【例4】 在坐标范围0 x2,-2y2内重新绘制正弦曲线，其程序为： x=linspace(0,2*pi,60); %生成含有60个数据元素的向量x y=sin(x); plot(x,y); axis (0 2*pi -2 2); %设定坐标轴范围,52,（四）加图例 给图形加图例命令为legend。该命令把图例放置在图形空白处，用户还可以通过鼠标移动图例，将其放到希望的位置。 格式:legend(图例说明,图例说明);,【例5】 为正弦、余弦曲线增加图例，其程序为： x=0:pi/100:2*pi; y1=sin(x); y2=cos(x); plot(x,y1,x,y2, -)

24、; legend(sin(x),cos(x);,53,（五）加网格线命令 若在图形中加网格线，用grid on。,阅读以下程序： x=-2:0.1:2; %产生横坐标x数组 y=x.3-3*x; %计算由y=x3-3x确定的纵坐标y数组 plot(x,y) %绘图 grid on %给图形加上网格线 axis equal %使x,y轴单位刻度相等,54,（一）subplot（m,n,p） 该命令将当前图形窗口分成mn个绘图区，即每行n个，共m行，区号按行优先编号，且选定第p个区为当前活动区。,二、 subplot 并列绘图命令,55,【例6】 在一个图形窗口中同时绘制正弦、余弦、正切、余切曲线

25、，程序为： x=linspace(0,2*pi,60); y=sin(x); z=cos(x); t=sin(x)./(cos(x)+eps); % eps为系统内部常数 ct=cos(x)./(sin(x)+eps); subplot(2,2,1); %分成22区域且指定1号为活动区 plot(x,y); title(sin(x); axis (0 2*pi -1 1); subplot(2,2,2); plot(x,z); title(cos(x); axis (0 2*pi -1 1); subplot(2,2,3); plot(x,t); title(tangent(x); axis

26、(0 2*pi -40 40); subplot(2,2,4); plot(x,ct); title(cotangent(x); axis (0 2*pi -40 40);,56,（二） figure 多图形窗口绘图命令 需要建立多个图形窗口，绘制并保持每一个窗口的图形，可以使用figure命令。 每执行一次figure命令，就创建一个新的图形窗口，该窗口自动为活动窗口，若需要还可以返回该窗口的识别号码，称该号码为句柄。句柄显示在图形窗口的标题栏中，即图形窗口标题。用户可通过句柄激活或关闭某图形窗口，而axis、xlabel、title等许多命令也只对活动窗口有效。,57,重新绘制上例4个图形

27、，程序变动后如下： x=linspace(0,2*pi,60); y=sin(x); z=cos(x); t=sin(x)./(cos(x)+eps); ct=cos(x)./(sin(x)+eps); H1=figure; %创建新窗口并返回句柄到变量H1 plot(x,y); %绘制图形并设置有关属性 title(sin(x); axis (0 2*pi -1 1); H2=figure; %创建第二个窗口并返回句柄到变量H2 plot(x,z); %绘制图形并设置有关属性 title(cos(x); axis (0 2*pi -1 1); H3=figure; %同上 plot(x,t)

28、; title(tangent(x); axis (0 2*pi -40 40); H4=figure; %同上 plot(x,ct); title(cotangent(x); axis (0 2*pi -40 40);,58,（三）hold图形保持命令 若在已存在图形窗口中用plot命令继续添加新的图形内容，可使用图形保持命令hold。发出命令hold on后，再执行plot命令，在保持原有图形或曲线的基础上，添加新绘制的图形。,59,阅读如下程序： x=linspace(0,2*pi,60); y=sin(x); z=cos(x); plot(x,y,b); %绘制正弦曲线 hold on

29、; %设置图形保持状态 plot(x,z,g); %保持正弦曲线同时绘制余弦曲线 axis (0 2*pi -1 1); legend(cos,sin); hold off %关闭图形保持,60,三、 fplot-函数f(x)绘图命令 fplot函数则可自适应地对函数进行采样，能更好地反应函数的变化规律。 fplot函数格式：fplot(fname，lims) 其中fname为函数名，以字符串形式出现， lims=a,b或a,b,c,d，为变量取值范围。 a,b为x的区间，c,d为y的区间。 例：fplot(sin(x),0 2*pi,-+) fplot(sin(x),cos(x),0 2*p

30、i,.) %同时绘制正弦、余弦曲线,61,为绘制f(x)=cos(tan(x)曲线，可先建立函数文件fct.m，其内容为： function y=fct(x) y=cos(tan(pi*x); 用fplot函数调用fct.m函数，其命令为： fplot(fct,0 1),62,四、 explot符号函数的绘图命令 ezplot函数格式：ezplot(fname，lims) 其中fname为函数名，以字符串形式出现， lims=a,b或a,b,c,d，为变量取值范围。 例：%绘制正弦函数从0到2pi区间上的图形ezplot(sin(x),0 2*pi) %绘制隐函数f(x,y)=0在a,b与c,

31、d区间上的图形 ezplot(4*x2+16*y2-3,-1 1 -1 1) %绘制参数方程x=sinx,y=cosx的图形 ezplot(sin(x),cos(x),0 2*pi),63,一、 对数坐标图形 （一）loglog(x,y) 双对数坐标 【例7】 绘制y=|1000sin(4x)|+1的双对数坐标图。程序为： x=0:0.1:2*pi; y=abs(1000*sin(4*x)+1; loglog(x,y); %双对数坐标绘图命令,2 特殊坐标图形,64,（二）单对数坐标 以X轴为对数重新绘制上述曲线，程序为： x=0:0.01:2*pi y=abs(1000*sin(4*x)+1

32、 semilogx(x,y); %单对数X轴绘图命令 同样，可以以Y轴为对数重新绘制上述曲线，程序为： x=0:0.01:2*pi y=abs(1000*sin(4*x)+1 semilogy(x,y); %单对数Y轴绘图命令,65,二、 极坐标图 函数polar(theta,rho)用来绘制极坐标图，theta为极坐标角度，rho为极坐标半径 【例8】 绘制sin(2*)*cos(2*)的极坐标图，程序为： theta=0:0.01:2*pi; rho=sin(2*theta).*cos(2*theta); polar(theta,rho); %绘制极坐标图命令 title(polar pl

33、ot);,66,一、阶梯图形 函数stairs(x,y)可以绘制阶梯图形，如下列程序段： x=-2.5:0.25:2.5; y=exp(-x.*x); stairs(x,y); %绘制阶梯图形命令 title(stairs plot);,3 其它图形函数,67,二、条形图形 函数bar(x,y)可以绘制条形图形，如下列程序段将绘制条形图形 x=-2.5:0.25:2.5; y=exp(-x.*x); bar(x,y); %绘制条形图命令,68,三、填充图形 fill(x,y,c)函数用来绘制并填充二维多边图形，x和y为二维多边形顶点坐标向量。字符 c 规定填充颜色，其取值前已叙述。 下述程序段

34、绘制一正方形并以黄色填充： x=0 1 1 0 0; %正方形顶点坐标向量 y=0 0 1 1 0; fill(x,y,y);%绘制并以黄色填充正方形图,69,再如： x=0:0.025:2*pi; y=sin(3*x); fill(x,y,0.5 0.3 0.4); %颜色向量 Matlab系统可用向量表示颜色，通常称其为颜色向量。基本颜色向量用r g b表示，即RGB颜色组合；以RGB为基本色，通过 r,g,b在01范围内的不同取值可以组合出各种颜色。,70,常用绘图命令： plot:用于数据点绘图。 fplot:用于函数绘图。 ezplot:用于符号函数绘图。可绘制隐函数和参数方程的图形

35、。 区别与差异： plot,fplot可对图形的线形，颜色作出控制，而ezplot则不能。 fplot可绘出比较精确的图形，而ezplot一般较适宜画不太精确的图形。,小结,71,plot 二维图形基本函数 fplot f(x)函数曲线绘制 ezplot 符号函数绘图 fill 填充二维多边图形 polar 极坐标图 bar 条形图 loglog 双对数坐标图 semilogx X轴为对数的坐标图 semilogy Y轴为对数的坐标图 stairs 阶梯形图,二维绘图函数小结,axis 设置坐标轴 figure 创建图形窗口 grid 放置坐标网格线 hold 保持当前图形窗口内容 subpl

36、ot 创建子图 title 放置图形标题 xlabel 放置X轴坐标标记 ylabel 放置Y轴坐标标记,72,阅读下面程序: %绘制摆线: hold on t=0:0.01:4*pi; for a=1:1:3 x=a*(t-sin(t); y=a*(1-cos(t); plot(x,y) end,73,h=3 2 1 0.5; %在曲线上取不同的点 a=(exp(h)-1)./h;%计算连接点M与与点P的各条割线的斜率 x=-1:0.1:3; %选定图形的自变量范围 plot(x,exp(x),r);%作函数图形 hold on; %在图形上继续作图 for i=1:4 plot(h(i),

37、exp(h(i),w) %在图上作出不同的点 plot(x,a(i)*x+1) %作割线的图 end axis square %把所有图形放在一个正方形框内 plot(x,x+1,g) %画出切线的图形,画出 在点P(0,1)处的切线及若干条割线,观察割线的变化趋势,理解导数的定义及几何意义.,74,一、 plot3函数 最基本的三维图形函数为plot3，它是将二维函数plot的有关功能扩展到三维空间，用来绘制三维图形。 函数格式：plot3(x1,y1,z1,c1,x2,y2,z2,c2,) 其中x1,y1,z1表示三维坐标向量，c1，c2表示线形或颜色。 函数功能：以向量x，y，z为坐标，

38、绘制三维曲线。,三维图形,75,【例9】 绘制三维螺旋曲线，其程序为： t=0:pi/50:10*pi; y1=sin(t);y2=cos(t); plot3(y1,y2,t); title(helix),text(0,0,0,origin); xlabel(sin(t),ylabel(cos(t),zlabel(t); grid;,76,二、mesh函数 mesh函数用于绘制三维网格图。在不需要绘制特别精细的三维曲面结构图时，可以通过绘制三维网格图来表示三维曲面。三维曲面的网格图最突出的优点是：它较好地解决了实验数据在三维空间的可视化问题。 函数格式：mesh(x,y,z,c) 其中x，y控

39、制X和Y轴坐标，矩阵z是由(x，y)求得Z轴坐标，(x,y,z)组成了三维空间的网格点；c用于控制网格点颜色。,【例10】 下列程序绘制三维网格曲面图 x=0:0.15:2*pi; y=0:0.15:2*pi; z=sin(y)*cos(x); %矩阵相乘 mesh(x,y,z);,77,三、surf函数 surf用于绘制三维曲面图，各线条之间的补面用颜色填充。surf函数和mesh函数的调用格式一致。 函数格式: surf (x,y,z) 其中x，y控制X和Y轴坐标，矩阵z是由x，y求得的曲面上Z轴坐标。,【例11】 下列程序绘制三维曲面图形 x=0:0.15:2*pi; y=0:0.15:

40、2*pi; z=sin(y)*cos(x); %矩阵相乘 surf(x,y,z); xlabel(x-axis),ylabel(y-axis),zlabel(z-label); title(3-D surf);,78,例绘制马鞍面的图形，并用平行截面法观察马鞍面的特点,x=-4:0.1:4; y=x; mx,my=meshgrid(x,y); mz=mx.2-my.2; ix=find(mx=2); px=2*ones(1,length(ix); py=my(ix); pz=mz(ix); subplot(1,2,1) hold on mesh(mx,my,mz) plot3(px,py,pz

41、,r*) subplot(1,2,2) plot3(px,py,pz),79,拟 合,2.拟合的基本原理,1. 拟合问题引例,80,拟 合 问 题 引 例 1,求600C时的电阻R。,设 R=at+b a,b为待定系数,81,拟 合 问 题 引 例 2,求血药浓度随时间的变化规律c(t).,作半对数坐标系(semilogy)下的图形,MATLAB(aa1),82,曲 线 拟 合 问 题 的 提 法,已知一组（二维）数据，即平面上 n个点（xi,yi) i=1,n, 寻求一个函数（曲线）y=f(x), 使 f(x) 在某种准则下与所有数据点最为接近，即曲线拟合得最好。,y=f(x),i 为点（x

42、i,yi) 与曲线 y=f(x) 的距离,83,拟合与插值的关系,函数插值与曲线拟合都是要根据一组数据构造一个函数作为近似，由于近似的要求不同，二者的数学方法上是完全不同的。,实例：下面数据是某次实验所得，希望得到X和 f之间的关系？,MATLAB(cn),问题：给定一批数据点，需确定满足特定要求的曲线或曲面,解决方案：,若不要求曲线（面）通过所有数据点，而是要求它反映对象整体的变化趋势，这就是数据拟合，又称曲线拟合或曲面拟合。,若要求所求曲线（面）通过所给所有数据点，就是插值问题；,84,最临近插值、线性插值、样条插值与曲线拟合结果：,85,曲线拟合问题最常用的解法线性最小二乘法的基本思路,

43、第一步:先选定一组函数 r1(x), r2(x), rm(x), mn, 令 f(x)=a1r1(x)+a2r2(x)+ +amrm(x) （1） 其中 a1,a2, am 为待定系数。,第二步: 确定a1,a2, am 的准则（最小二乘准则）： 使n个点（xi,yi) 与曲线 y=f(x) 的距离i 的平方和最小 。,记,问题归结为，求 a1,a2, am 使 J(a1,a2, am) 最小。,86,线性最小二乘法的求解：预备知识,超定方程组：方程个数大于未知量个数的方程组,超定方程一般是不存在解的矛盾方程组。,如果有向量a使得 达到最小， 则称a为上述超定方程的最小二乘解。,87,线性最小

44、二乘法的求解,定理：当RTR可逆时，超定方程组（3）存在最小二乘解，且即为方程组 RTRa=RTy 的解：a=(RTR)-1RTy,所以，曲线拟合的最小二乘法要解决的问题，实际上就是求以下超定方程组的最小二乘解的问题。,88,线性最小二乘拟合 f(x)=a1r1(x)+ +amrm(x)中函数r1(x), rm(x)的选取,1. 通过机理分析建立数学模型来确定 f(x)；,2. 将数据 (xi,yi) i=1, n 作图，通过直观判断确定 f(x)：,89,用MATLAB解拟合问题,1、线性最小二乘拟合,2、非线性最小二乘拟合,90,用MATLAB作线性最小二乘拟合,1. 作多项式f(x)=a

45、1xm+ +amx+am+1拟合,可利用已有程序:,a=polyfit(x,y,m),2. 对超定方程组,3.多项式在x处的值y可用以下命令计算： y=polyval（a，x）,91,例 对下面一组数据作二次多项式拟合,92,1）输入以下命令： x=0:0.1:1; y=-0.447 1.978 3.28 6.16 7.08 7.34 7.66 9.56 9.48 9.30 11.2; R=(x.2) x ones(11,1)； A=Ry,MATLAB(zxec1),解法1用解超定方程的方法,2）计算结果： = -9.8108 20.1293 -0.0317,93,1）输入以下命令： x=0:

46、0.1:1; y=-0.447 1.978 3.28 6.16 7.08 7.34 7.66 9.56 9.48 9.30 11.2; A=polyfit(x,y,2) z=polyval(A,x); plot(x,y,k+,x,z,r) %作出数据点和拟合曲线的图形,2）计算结果： = -9.8108 20.1293 -0.0317,解法2用多项式拟合的命令,MATLAB(zxec2),94,1. lsqcurvefit 已知数据点： xdata=（xdata1，xdata2，xdatan）， ydata=（ydata1，ydata2，ydatan）,用MATLAB作非线性最小二乘拟合,Ma

47、tlab的提供了两个求非线性最小二乘拟合的函数：lsqcurvefit和lsqnonlin。两个命令都要先建立M-文件fun.m，在其中定义函数f(x)，但两者定义f(x)的方式是不同的,可参考例题.,lsqcurvefit用以求含参量x（向量）的向量值函数 F(x,xdata)=（F（x，xdata1），F（x，xdatan）T 中的参变量x(向量),使得,95,输入格式为: (1) x = lsqcurvefit (fun,x0,xdata,ydata); (2) x =lsqcurvefit (fun,x0,xdata,ydata,options); (3) x = lsqcurvefi

48、t (fun,x0,xdata,ydata,options,grad); (4) x, options = lsqcurvefit (fun,x0,xdata,ydata,); (5) x, options,funval = lsqcurvefit (fun,x0,xdata,ydata,); (6) x, options,funval, Jacob = lsqcurvefit (fun,x0,xdata,ydata,);,说明：x = lsqcurvefit (fun,x0,xdata,ydata,options);,96,lsqnonlin用以求含参量x（向量）的向量值函数 f(x)=(f

49、1(x),f2(x),fn(x)T 中的参量x，使得 最小。 其中 fi（x）=f（x，xdatai，ydatai） =F(x,xdatai)-ydatai,2. lsqnonlin,已知数据点： xdata=（xdata1，xdata2，xdatan） ydata=（ydata1，ydata2，ydatan）,97,输入格式为： 1） x=lsqnonlin（fun，x0）； 2） x= lsqnonlin （fun，x0，options）； 3） x= lsqnonlin （fun，x0，options，grad）； 4） x，options= lsqnonlin （fun，x0，）； 5

50、） x，options，funval= lsqnonlin （fun，x0，）；,说明：x= lsqnonlin （fun，x0，options）；,98,例2 用下面一组数据拟合 中的参数a，b，k,该问题即解最优化问题：,99,MATLAB(fzxec1),1）编写M-文件 curvefun1.m function f=curvefun1(x,tdata) f=x(1)+x(2)*exp(-0.02*x(3)*tdata) %其中 x(1)=a; x(2)=b；x(3)=k;,2）输入命令 tdata=100:100:1000 cdata=1e-03*4.54,4.99,5.35,5.65

51、,5.90,6.10,6.26,6.39, 6.50,6.59; x0=0.2,0.05,0.05; x=lsqcurvefit (curvefun1,x0,tdata,cdata) f= curvefun1(x,tdata),F(x，tdata)= ，x=(a，b，k),解法1. 用命令lsqcurvefit,100,3）运算结果为： f =0.0043 0.0051 0.0056 0.0059 0.0061 0.0062 0.0062 0.0063 0.0063 0.0063 x = 0.0063 -0.0034 0.2542,4）结论：a=0.0063, b=-0.0034, k=0.2

52、542,101,MATLAB(fzxec2),解法 2 用命令lsqnonlin f(x)=F(x,tdata,ctada)= x=（a，b，k）,1）编写M-文件 curvefun2.m function f=curvefun2(x) tdata=100:100:1000; cdata=1e-03*4.54,4.99,5.35,5.65,5.90, 6.10,6.26,6.39,6.50,6.59; f=x(1)+x(2)*exp(-0.02*x(3)*tdata)- cdata,2）输入命令: x0=0.2,0.05,0.05; x=lsqnonlin(curvefun2,x0) f= c

53、urvefun2(x),函数curvefun2的自变量是x，cdata和tdata是已知参数，故应将cdata tdata的值写在curvefun2.m中,102,3）运算结果为 f =1.0e-003 *(0.2322 -0.1243 -0.2495 -0.2413 -0.1668 -0.0724 0.0241 0.1159 0.2030 0.2792 x =0.0063 -0.0034 0.2542,可以看出,两个命令的计算结果是相同的.,4）结论：即拟合得a=0.0063 b=-0.0034 k=0.2542,103,MATLAB解应用问题实例,1、电阻问题,2、给药方案问题,*3、水塔

54、流量估计问题,104,MATLAB(dianzu1),电阻问题,得到 a1=3.3940, a2=702.4918,方法2.直接用,结果相同。,MATLAB(dianzu2),105,一室模型：将整个机体看作一个房室，称中心室，室内血药浓度是均匀的。快速静脉注射后，浓度立即上升；然后迅速下降。当浓度太低时，达不到预期的治疗效果；当浓度太高，又可能导致药物中毒或副作用太强。临床上，每种药物有一个最小有效浓度c1和一个最大有效浓度c2。设计给药方案时，要使血药浓度 保持在c1c2之间。本题设c1=10，c2=25(ug/ml).,一种新药用于临床之前，必须设计给药方案.,药物进入机体后血液输送到全

55、身，在这个过程中不断地被吸收、分布、代谢，最终排出体外，药物在血液中的浓度，即单位体积血液中的药物含量，称为血药浓度。,106,要设计给药方案,必须知道给药后血药浓度随时间变化的规律。从实验和理论两方面着手：,107,给药方案,1. 在快速静脉注射的给药方式下，研究血药浓度（单位体积血液中的药物含量）的变化规律。,t,问题,2. 给定药物的最小有效浓度和最大治疗浓度，设计给药方案：每次注射剂量多大；间隔时间多长。,分析,理论：用一室模型研究血药浓度变化规律,实验：对血药浓度数据作拟合，符合负指数变化规律,108,3.血液容积v, t=0注射剂量d, 血药浓度立即为d/v.,2.药物排除速率与血

56、药浓度成正比，比例系数 k(0),模型假设,1. 机体看作一个房室，室内血药浓度均匀一室模型,模型建立,在此，d=300mg，t及c（t）在某些点处的值见前表，需经拟合求出参数k、v,109,用线性最小二乘拟合c(t),MATLAB(lihe1),计算结果：,用非线性最小二乘拟合c(t),110,给药方案 设计,设每次注射剂量D, 间隔时间,血药浓度c(t) 应c1 c(t) c2,初次剂量D0 应加大,给药方案记为：,2、,1、,计算结果：,给药方案：,c1=10,c2=25 k=0.2347 v=15.02,111,故可制定给药方案：,即: 首次注射375mg， 其余每次注射225mg，

57、注射的间隔时间为4小时。,112,估计水塔的流量,2、解题思路,3、算法设计与编程,1、问题,113,某居民区有一供居民用水的园柱形水塔，一般可以通过测量其水位来估计水的流量，但面临的困难是，当水塔水位下降到设定的最低水位时，水泵自动启动向水塔供水，到设定的最高水位时停止供水，这段时间无法测量水塔的水位和水泵的供水量通常水泵每天供水一两次，每次约两小时. 水塔是一个高12.2米，直径17.4米的正园柱按照设计，水塔水位降至约8.2米时，水泵自动启动，水位升到约10.8米时水泵停止工作 表1 是某一天的水位测量记录，试估计任何时刻（包括水泵正供水时）从水塔流出的水流量，及一天的总用水量,114,

58、115,流量估计的解题思路,拟合水位时间函数,确定流量时间函数,估计一天总用水量,116,拟合水位时间函数 测量记录看，一天有两个供水时段（以下称第1供水时段和第2供水时段），和3个水泵不工作时段（以下称第1时段t=0到t=8.97，第2次时段t=10.95到t=20.84和第3时段t=23以后）对第1、2时段的测量数据直接分别作多项式拟合，得到水位函数为使拟合曲线比较光滑，多项式次数不要太高，一般在36由于第3时段只有3个测量记录，无法对这一时段的水位作出较好的拟合,117,2、确定流量时间函数 对于第1、2时段只需将水位函数求导数即可，对于两个供水时段的流量，则用供水时段前后（水泵不工作时段）的流量拟合得到，并且将拟合得到的第2供水时段流量外推，将第3时段流量包含在第2供水时段内,118,3、一天总用水