MATLAB在数学建模中的应用.ppt

《MATLAB在数学建模中的应用.ppt》由会员分享，可在线阅读，更多相关《MATLAB在数学建模中的应用.ppt（131页珍藏版）》请在三一办公上搜索。

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

2、能.在国际学术界，MATLAB已经被确认为准确、可靠的科学计算标准软件.在设计研究单位和工业部门，MATLAB被认作进行高效研究、开发的首选软件工具.,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版,MATLAB的功能,MATLAB产品组是从支持概念设计、算法开发、建模仿真，到实时实现的集成环境，可

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

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

5、在微积分中的应用,1、求函数值,例1 在命令窗口中键入表达式并求 时的函数值。,x=2,y=4z=x2+exp(x+y)-y*log(x)-3,x=2y=4z=401.6562,命令窗口显示结果：,例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;endy,y=e(1,100)ans=y y=2.7183,调用函数M文件,MATLAB在微积分中的应用,2、求极限,例3 求极限,syms n;limit(sqrt(n+sqrt(n)-sqrt(n),n,inf),a

6、ns=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.,定义符号变量,MATLAB在微积分中的应用,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为符号变量,求,Dif

7、ference：差分 Differential：微分的,syms x;y=log(1+x);a=diff(y,x,2)a=-1/(1+x)2 x=1;eval(a)ans=-0.2500,求,求,将符号表达式转换成数值表达式,例6 设,，求,syms x y;z=exp(2*x)*(x+y2+2*y);a=diff(z,x)b=diff(z,y)c=diff(z,x,2)d=diff(z,y,2)e=diff(a,y),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*

8、exp(2*x)e=2*exp(2*x)*(2*y+2),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-005,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.0008FVAL=-2.5000,MATLAB在微积分中的应用,5、求积分,例9 求不定积分

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(1);y=x.2.*log(x);trapz(x,y)ans=4.5137,例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

10、)ans=0.8556,变步长数值积分,梯形法数值积分,MATLAB在微积分中的应用,5、求积分,例11 求二重积分,syms x y;f=y2/x2;int(int(f,x,1/2,2),y,1,2)ans=7/2,符号积分,f=(y.2)./(x.2);dblquad(f,1/2,2,1,2)ans=3.5000,数值计算,MATLAB在微积分中的应用,6、解微分方程,例12 计算初值问题：,dsolve(Dy=x+y,y(0)=1,x),ans=-x-1+2*exp(x),一定要大写,MATLAB在微积分中的应用,7、级数问题,例13 求函数 的泰勒展开式，并计算该函数在x=3.42时的

11、近似值。,syms x;taylor(sin(x)/x,x,10),ans=1-1/6*x2+1/120*x4-1/5040*x6+1/362880*x8,x=3.42;eval(ans)ans=-0.0753,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,MATLAB在线性代数中的应用,1、矩阵的基本运算,例1 已知,inv(a)ans=0.1582-0.1013 0.0633-0.0886-0.0633 0.1646 0.0949 0.1392 0.0380,MATLAB在

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

13、 0 0 1,将矩阵A化为最简阶梯形,R（A）=4=n；所以方程组只有零解。,RREF Reduced row echelon form,2、解线性方程组,求齐次方程组的基础解系,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,所以方程组的一般解为,3、将矩阵对角化,a=-1 2 0;-2 3 0;3 0 2;v,d=eig(a)v=0 379/1257 379/1257 0 379/1257 379/

14、1257 1-379/419-379/419 d=2 0 0 0 1 0 0 0 1,A的特征值为2，1，1,4、用正交变换化二次型为标准形,a=1 1 1 11 1 1 11 1 1 11 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.5000t=-0.0000 0 0 0 0-0.0000 0 0 0 0 0 0 0 0 0 4.0000,a=1 1 1 1；1 1 1 1；1 1 1 1

15、；1 1 1 1;format ratu 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 integers.,4、用正交变换化二次型为标准形,结论：作正交变换,则有,上机实验题一、基础型实验1、计算下列极限,2、计算下列导数

16、(1)(2)(3)(4),实验练习,一.输入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,-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

17、=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.,四.将下列矩阵化为阶梯矩阵,1.A=1,-2,0;-1,1,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.,五.求下列矩阵的

18、秩 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.,附录：MATLAB软件中部分常用函数表,作为一个功能强大的工具软件，Matlab具有很强的图形处理功能，提供了大量的二维、三维图形函数。由于系统采用面向对象的技术和丰富的矩阵运算，所以在图形处理方面即方便又高效。,一、plot数据点绘图命令命令格式：plot(x,y)其中x和y为坐标向量命令功能：以向量x、y为轴，绘制曲线。【例1】在区

19、间0X2内，绘制正弦曲线Y=sin(x)，其程序为：x=0:pi/100:2*pi;y=sin(x);plot(x,y),1.二维图形,【例2】同时绘制正、余弦两条曲线y1=sin(x)和y2=cos(x)，其程序为：x=0:pi/100:2*pi;y1=sin(x);y2=cos(x);plot(x,y1,x,y2)plot函数还可以为plot(x,y1,x,y2，x,y3，)形式，其功能是以公共向量x为X轴，分别以y1，y2，y3，为Y轴，在同一幅图内绘制出多条曲线。,（一）线型与颜色格式：plot(x,y1,cs,.)其中c表示颜色，s表示线型。,【例3】用不同线型和颜色重新绘制例2图形

20、，其程序为：x=0:pi/100:2*pi;y1=sin(x);y2=cos(x);plot(x,y1,go,x,y2,b-.)其中参数go和b-.表示图形的颜色和线型。g表示绿色，o表示图形线型为圆圈；b表示蓝色，-.表示图形线型为点划线。,绘图基本线型和颜色,（二）图形标记在绘制图形的同时，可以对图形加上一些说明，如图形名称、图形某一部分的含义、坐标说明等，将这些操作称为添加图形标记。title(加图形标题);xlabel(加X轴标记);ylabel(加Y轴标记);text(X,Y,添加文本);,（三）设定坐标轴用户若对坐标系统不满意，可利用axis命令对其重新设定。axis(xmin x

21、max ymin ymax)设定最大和最小值axis（auto）将坐标系统返回到自动缺省状态axis（square）将当前图形设置为方形axis（equal）两个坐标因子设成相等axis（off）关闭坐标系统axis（on）显示坐标系统,【例4】在坐标范围0 x2,-2y2内重新绘制正弦曲线，其程序为：x=linspace(0,2*pi,60);%生成含有60个数据元素的向量xy=sin(x);plot(x,y);axis(0 2*pi-2 2);%设定坐标轴范围,（四）加图例给图形加图例命令为legend。该命令把图例放置在图形空白处，用户还可以通过鼠标移动图例，将其放到希望的位置。格式:l

22、egend(图例说明,图例说明);,【例5】为正弦、余弦曲线增加图例，其程序为：x=0:pi/100:2*pi;y1=sin(x);y2=cos(x);plot(x,y1,x,y2,-);legend(sin(x),cos(x);,（五）加网格线命令若在图形中加网格线，用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轴单位刻度相等,（一）subplot（m,n,p）该命令将当前图形窗口分成mn个绘图区，即每行n个，共m行，区

23、号按行优先编号，且选定第p个区为当前活动区。,二、subplot 并列绘图命令,【例6】在一个图形窗口中同时绘制正弦、余弦、正切、余切曲线，程序为：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);su

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

25、le等许多命令也只对活动窗口有效。,重新绘制上例4个图形，程序变动后如下：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;%创建新窗口并返回句柄到变量H1plot(x,y);%绘制图形并设置有关属性title(sin(x);axis(0 2*pi-1 1);H2=figure;%创建第二个窗口并返回句柄到变量H2plot(x,z);%绘制图形并设置有关属性title(cos(x);axis(0 2*pi-1 1);H3=figure;%同上plot(

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

27、);%保持正弦曲线同时绘制余弦曲线axis(0 2*pi-1 1);legend(cos,sin);hold off%关闭图形保持,三、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*pi,.)%同时绘制正弦、余弦曲线,为绘制f(x)=cos(tan(x)曲线，可先

28、建立函数文件fct.m，其内容为：function y=fct(x)y=cos(tan(pi*x);用fplot函数调用fct.m函数，其命令为：fplot(fct,0 1),四、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,d区间上的图形ezplot(4*x2+16*y2-3,-1 1-1 1)%绘制参数方程x=sinx,y=cosx的

29、图形ezplot(sin(x),cos(x),0 2*pi),一、对数坐标图形（一）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 特殊坐标图形,（二）单对数坐标以X轴为对数重新绘制上述曲线，程序为：x=0:0.01:2*piy=abs(1000*sin(4*x)+1semilogx(x,y);%单对数X轴绘图命令同样，可以以Y轴为对数重新绘制上述曲线，程序为：x=0:0.01:2*piy=abs(1000*sin(4*x

30、)+1semilogy(x,y);%单对数Y轴绘图命令,二、极坐标图函数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 plot);,一、阶梯图形函数stairs(x,y)可以绘制阶梯图形，如下列程序段：x=-2.5:0.25:2.5;y=exp(-x.*x);stairs(x,y);%绘制阶梯图形命令title(st

31、airs plot);,3 其它图形函数,二、条形图形函数bar(x,y)可以绘制条形图形，如下列程序段将绘制条形图形x=-2.5:0.25:2.5;y=exp(-x.*x);bar(x,y);%绘制条形图命令,三、填充图形fill(x,y,c)函数用来绘制并填充二维多边图形，x和y为二维多边形顶点坐标向量。字符 c 规定填充颜色，其取值前已叙述。下述程序段绘制一正方形并以黄色填充：x=0 1 1 0 0;%正方形顶点坐标向量y=0 0 1 1 0;fill(x,y,y);%绘制并以黄色填充正方形图,再如：x=0:0.025:2*pi;y=sin(3*x);fill(x,y,0.5 0.3 0

32、.4);%颜色向量Matlab系统可用向量表示颜色，通常称其为颜色向量。基本颜色向量用r g b表示，即RGB颜色组合；以RGB为基本色，通过 r,g,b在01范围内的不同取值可以组合出各种颜色。,常用绘图命令：plot:用于数据点绘图。fplot:用于函数绘图。ezplot:用于符号函数绘图。可绘制隐函数和参数方程的图形。区别与差异：plot,fplot可对图形的线形，颜色作出控制，而ezplot则不能。fplot可绘出比较精确的图形，而ezplot一般较适宜画不太精确的图形。,小结,plot 二维图形基本函数fplot f(x)函数曲线绘制ezplot 符号函数绘图fill 填充二维多边图

33、形polar 极坐标图bar 条形图loglog 双对数坐标图semilogx X轴为对数的坐标图semilogy Y轴为对数的坐标图stairs 阶梯形图,二维绘图函数小结,axis 设置坐标轴figure 创建图形窗口grid 放置坐标网格线hold 保持当前图形窗口内容subplot 创建子图title 放置图形标题xlabel 放置X轴坐标标记ylabel 放置Y轴坐标标记,阅读下面程序:%绘制摆线:hold ont=0:0.01:4*pi;for a=1:1:3 x=a*(t-sin(t);y=a*(1-cos(t);plot(x,y)end,h=3 2 1 0.5;%在曲线上取不同

34、的点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),exp(h(i),w)%在图上作出不同的点 plot(x,a(i)*x+1)%作割线的图endaxis square%把所有图形放在一个正方形框内plot(x,x+1,g)%画出切线的图形,画出 在点P(0,1)处的切线及若干条割线,观察割线的变化趋势,理解导数的定义及几何意义.,一、plot3函数最基本的三维图形函数为plot3，它是将二维函数plot的有关

35、功能扩展到三维空间，用来绘制三维图形。函数格式：plot3(x1,y1,z1,c1,x2,y2,z2,c2,)其中x1,y1,z1表示三维坐标向量，c1，c2表示线形或颜色。函数功能：以向量x，y，z为坐标，绘制三维曲线。,三维图形,【例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;,二、mesh函数mesh函数用于绘制三维网格图。在不需要绘制特别精细的三维曲

36、面结构图时，可以通过绘制三维网格图来表示三维曲面。三维曲面的网格图最突出的优点是：它较好地解决了实验数据在三维空间的可视化问题。函数格式：mesh(x,y,z,c)其中x，y控制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);,三、surf函数surf用于绘制三维曲面图，各线条之间的补面用颜色填充。surf函数和mesh函数的调用格式一致。函数格式:surf(x,y,z)其中x

37、，y控制X和Y轴坐标，矩阵z是由x，y求得的曲面上Z轴坐标。,【例11】下列程序绘制三维曲面图形x=0:0.15:2*pi;y=0:0.15:2*pi;z=sin(y)*cos(x);%矩阵相乘surf(x,y,z);xlabel(x-axis),ylabel(y-axis),zlabel(z-label);title(3-D surf);,例绘制马鞍面的图形，并用平行截面法观察马鞍面的特点,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(i

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

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

40、问题最常用的解法线性最小二乘法的基本思路,第一步:先选定一组函数 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)最小。,线性最小二乘法的求解：预备知识,超定方程组：方程个数大于未知量个数的方程组,超定方程一般是不存在解的矛盾方程组。,如果有向量a使得 达到最小，则称a为上述超定方程的最小二乘解。,线性最小二乘法的

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

42、,a=polyfit(x,y,m),2.对超定方程组,3.多项式在x处的值y可用以下命令计算：y=polyval（a，x）,例 对下面一组数据作二次多项式拟合,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,1）输入以下命令：x=0:0.1:1;y=-0.447 1.978 3.28 6.16 7.08 7.34 7.66

43、 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),1.lsqcurvefit已知数据点：xdata=（xdata1，xdata2，xdatan），ydata=（ydata1，ydata2，ydatan）,用MATLAB作非线性最小二乘拟合,Matlab的提供了两个求非线性最小二乘拟合的函数：lsqcurvefit和lsqnonlin。两个命令都要先建立M-文件f

44、un.m，在其中定义函数f(x)，但两者定义f(x)的方式是不同的,可参考例题.,lsqcurvefit用以求含参量x（向量）的向量值函数F(x,xdata)=（F（x，xdata1），F（x，xdatan）T中的参变量x(向量),使得,输入格式为:(1)x=lsqcurvefit(fun,x0,xdata,ydata);(2)x=lsqcurvefit(fun,x0,xdata,ydata,options);(3)x=lsqcurvefit(fun,x0,xdata,ydata,options,grad);(4)x,options=lsqcurvefit(fun,x0,xdata,ydata

45、,);(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);,lsqnonlin用以求含参量x（向量）的向量值函数 f(x)=(f1(x),f2(x),fn(x)T 中的参量x，使得 最小。其中 fi（x）=f（x，xdatai，ydatai）=F(x,xdatai)-ydatai,2.lsqnonlin,已知数据点：xdata=（xd

46、ata1，xdata2，xdatan）ydata=（ydata1，ydata2，ydatan）,输入格式为：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）x，options，funval=lsqnonlin（fun，x0，）；,说明：x=lsqnonlin（fun，x0，options）；,例2 用下面一组数据拟合 中的参数a，b，k,该问题即解最优化问题：,MATLAB(fzxec1),1）编写M-文

47、件 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:1000cdata=1e-03*4.54,4.99,5.35,5.65,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.用命令lsqcur

48、vefit,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.2542,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.6

49、5,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=curvefun2(x),函数curvefun2的自变量是x，cdata和tdata是已知参数，故应将cdata tdata的值写在curvefun2.m中,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

50、-0.0034 0.2542,可以看出,两个命令的计算结果是相同的.,4）结论：即拟合得a=0.0063 b=-0.0034 k=0.2542,MATLAB解应用问题实例,1、电阻问题,2、给药方案问题,*3、水塔流量估计问题,MATLAB(dianzu1),电阻问题,得到 a1=3.3940,a2=702.4918,方法2.直接用,结果相同。,MATLAB(dianzu2),一室模型：将整个机体看作一个房室，称中心室，室内血药浓度是均匀的。快速静脉注射后，浓度立即上升；然后迅速下降。当浓度太低时，达不到预期的治疗效果；当浓度太高，又可能导致药物中毒或副作用太强。临床上，每种药物有一个最小有效