MATLAB在数学建模中的应用学习教案

1、会计学1MATLAB在数学在数学(shxu)建模中的应用建模中的应用第一页，共132页。第1页/共131页第二页，共132页。问题(wnt)的分析修正(xizhng)模型粗假设修正算法结果分析讨论推广修正假设粗模型粗算法发现问题发现规律模型验证第2页/共131页第三页，共132页。 Matlab软件(run jin)简介 数学建模Matlab算法第3页/共131页第四页，共132页。第4页/共131页第五页，共132页。第5页/共131页第六页，共132页。MATLABMATLAB产品组是从支持概念设计、算法开发、建模仿真，产品组是从支持概念设计、算法开发、建模仿真，到实时实现的集成环境，可用

2、来进行：到实时实现的集成环境，可用来进行：数据分析数据分析数值数值(shz)(shz)与符号计算与符号计算工程与科学绘图工程与科学绘图控制系统设计控制系统设计数字图像信号处理数字图像信号处理建模、仿真、原型开发建模、仿真、原型开发财务工程、应用开发、图形用户界面设计财务工程、应用开发、图形用户界面设计第6页/共131页第七页，共132页。编程效率高，允许编程效率高，允许(ynx)(ynx)用数学的语言来编写程序用数学的语言来编写程序用户使用方便，把程序的编辑、编译、连接和执行融为一体用户使用方便，把程序的编辑、编译、连接和执行融为一体高效方便的矩阵和数组运算高效方便的矩阵和数组运算语句简单，内

3、涵丰富语句简单，内涵丰富扩充能力强，交互性，开放性扩充能力强，交互性，开放性方便的绘图功能方便的绘图功能该软件由该软件由c c语言编写，移植性好语言编写，移植性好第7页/共131页第八页，共132页。第8页/共131页第九页，共132页。第9页/共131页第十页，共132页。第10页/共131页第十一页，共132页。菜单项；菜单项；工具栏；工具栏；【Command WindowCommand Window】命令窗口；】命令窗口；【Launch PadLaunch Pad】分类帮助窗口；】分类帮助窗口；【WorkspaceWorkspace】工作】工作(gngzu)(gngzu)区窗口；区窗口；

4、【Command HistoryCommand History】指令历史记录窗口；】指令历史记录窗口；【Current DirectoryCurrent Directory】当前目录选择窗口；】当前目录选择窗口；第11页/共131页第十二页，共132页。l MATLAB操作操作(cozu)窗口窗口双击桌面快捷键，启动软件。双击桌面快捷键，启动软件。接受命令的窗口接受命令的窗口(chungku)第12页/共131页第十三页，共132页。MATLAB在微积分中的应用在微积分中的应用(yngyng) 1、求函数值、求函数值 例例1 在命令窗口中键入表达式在命令窗口中键入表达式并求并求 时的函数值。时

5、的函数值。2ln3,x yzxeyx2,4xy x=2,y=4z=x2+exp(x+y)-y*log(x)-3x = 2y = 4z = 401.6562命令命令(mng lng)窗口显示结果：窗口显示结果： 第13页/共131页第十四页，共132页。例例2 用循环语句编写用循环语句编写M文件文件(wnjin)计算计算ex的值，其中的值，其中x,n为输入为输入变量，变量，ex的近似表达式为的近似表达式为2312!3!nxxxxexn 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.

6、7183调用函数调用函数M文件文件(wnjin)第14页/共131页第十五页，共132页。MATLAB在微积分中的应用在微积分中的应用(yngyng) 2、求极限、求极限(jxin) 例例3 求极限求极限(jxin) lim()nnnn syms n; limit(sqrt(n+sqrt(n)-sqrt(n),n,inf) ans = 1/2LIMIT 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,

7、a,left) specify the direction of a one-sided limit.定义符号变量定义符号变量第15页/共131页第十六页，共132页。MATLAB在微积分中的应用在微积分中的应用(yngyng) 3、求导数、求导数(do sh) 例例4 设设 1010lnxyxx，求，求 y 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为符号为符号(fho)变量变量 求求 dydxDifference：差分 Differential：微分的 第16页/

8、共131页第十七页，共132页。例例5 设设 ln(1),yx求求 212xd ydx syms x; y=log(1+x); a=diff(y,x,2) a = -1/(1+x)2 x=1;eval(a)ans = -0.2500求求 22d ydx求求 212xd ydx将符号将符号(fho)表达式表达式转换成数值表达式转换成数值表达式第17页/共131页第十八页，共132页。例例6 设设 222xzexyy，求，求 22222,zzzzzxyxyx y 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=

9、diff(z,y,2)e=diff(a,y) zaxzby22zcx22zdy2azeyx y 第18页/共131页第十九页，共132页。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)222222xxzaexyyex222xzbeyy22222424xxzcexyyex2222xzdey22222xzeeyx y 第19页/共131页第二十页，共132页。MATLAB在微积分中的应用在微积分中

10、的应用(yngyng) 4、求极值、求极值(j zh)和零点和零点 例例7 已知已知 5432( )323f xxxxx，求，求 （1）函数的零点；（）函数的零点；（2）函数在）函数在-1，2上的最小值上的最小值 fzero(3*x5-x4+2*x3+x2+3,0)ans = -0.8952 起始起始(q sh)点点 函数函数 命令函数命令函数 fminbnd(3*x5-x4+2*x3+x2+3,-1,2)ans = -1.1791e-005第20页/共131页第二十一页，共132页。MATLAB在微积分中的应用在微积分中的应用(yngyng) 4、求极值、求极值(j zh)和零点和零点 ，求

11、，求 例例8 已知已知 222( , , )2.5sinf x y zxyxy z函数在点（函数在点（1，-1，0）附近的最小值）附近的最小值 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第21页/共131页第二十二页，共132页。MATLAB在微积分中的应用在微积分中的应用(yngyng) 5、求积分、求积分(jfn) 例例9 求不定积分求不定积分(b dn j fn) cos2 cos3xxdx int(cos(2*x)*cos(3*x

12、) ans =1/2*sin(x)+1/10*sin(5*x)例例10 求定积分求定积分 21lnexxdxIntegrate：积分 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第22页/共131页第二十三页，共132页。例例10 求定积分求定积分(jfn) 2120 xedx int(exp(-x2/2),0,1) ans = 1/2*erf(1/2*2(1/2)*2(1/2)*pi(1/2)202( )xterf xedt22 2022tansed

13、t 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变步长变步长数值积分数值积分 梯形梯形(txng)法数值积分法数值积分 第23页/共131页第二十四页，共132页。MATLAB在微积分中的应用在微积分中的应用(yngyng) 5、求积分、求积分(jfn) 例例11 求二重积分求二重积分 221,2,122ydxdyxyx syms x y; f=y2/x2; int(int(f,x,1/2,2),y,1,2) ans =7/2符号符号(fho)积分积分 f=(y.2

14、)./(x.2); dblquad(f,1/2,2,1,2)ans = 3.5000数值计算数值计算 第24页/共131页第二十五页，共132页。MATLAB在微积分中的应用在微积分中的应用(yngyng) 6、解微分方程、解微分方程(wi fn fn chn) 例12 计算(j sun)初值问题：1)0(yxydxdy dsolve(Dy=x+y,y(0)=1,x)ans =-x-1+2*exp(x)一定要大写一定要大写 第25页/共131页第二十六页，共132页。MATLAB在微积分中的应用在微积分中的应用(yngyng) 7、级数、级数(j sh)问题问题 例例13 求函数求函数 的泰勒

15、展开式，并计算该的泰勒展开式，并计算该函数在函数在x=3.42时的近似值。时的近似值。sin( )xf xx 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第26页/共131页第二十七页，共132页。MATLAB在线性代数在线性代数(xin xn di sh)中的应用中的应用 1、矩阵、矩阵(j zhn)的基本运算的基本运算 例例1 已知已知 422134305 ,203153211AB a=4 -2 2;-3 0 5;1 5 3;

16、b=1 3 4;-2 0 -3;2 -1 1; a*b12 10 24 7 -14 -7-3 0 -8ans =AB 第27页/共131页第二十八页，共132页。MATLAB在线性代数在线性代数(xin xn di sh)中的应用中的应用 1、矩阵、矩阵(j zhn)的基本运算的基本运算 例例1 已知已知 422134305 ,203153211AB inv(a)ans = 0.1582 -0.1013 0.0633 -0.0886 -0.0633 0.1646 0.0949 0.1392 0.03801A第28页/共131页第二十九页，共132页。MATLAB在线性代数在线性代数(xin x

17、n di sh)中的应用中的应用 1、矩阵的基本、矩阵的基本(jbn)运算运算 例例1 已知已知 422134305 ,203153211AB ( )R A rank(a)ans = 3第29页/共131页第三十页，共132页。MATLAB在线性代数在线性代数(xin xn di sh)中的应用中的应用 1、矩阵的基本、矩阵的基本(jbn)运算运算 例例1 已知已知 422134305 ,203153211AB 1AB a/bans = 0 0 2.0000 -2.7143 -8.0000 -8.1429 2.4286 3.0000 2.2857第30页/共131页第三十一页，共132页。MA

18、TLAB在线性代数在线性代数(xin xn di sh)中的应用中的应用 1、矩阵、矩阵(j zhn)的基本运算的基本运算 例例1 已知已知 422134305 ,203153211AB 1A B abans = 0.4873 0.4114 1.0000 0.3671 -0.4304 0 -0.1076 0.2468 0第31页/共131页第三十二页，共132页。2、解线性方程组、解线性方程组 123412341234123442020372031260 xxxxxxxxxxxxxxxx a=1 -1 4 -2;1 -1 -1 2;3 1 7 -2;1 -3 -12 6; rref(a)ans

19、 =1 0 0 00 1 0 00 0 1 00 0 0 1将矩阵将矩阵(j zhn)A化为最简阶梯形化为最简阶梯形R（A）=4=n；所以所以(suy)方程组只有零解。方程组只有零解。RREF Reduced row echelon form第32页/共131页第三十三页，共132页。2、解线性方程组、解线性方程组 12312312312323424538213496xxxxxxxxxxxx 第33页/共131页第三十四页，共132页。求齐次方程组求齐次方程组的基础的基础(jch)解解系系 a=2 3 1;1 -2 4;3 8 -2;4 -1 9; b=4;-5;13;-6; c=null(a

20、,r)c = -2 1 1 求非齐次方程组求非齐次方程组的一个的一个(y )特解特解 l u=lu(a); x0=u(lb)x0 = -3124/135 3529/270 2989/270 所以所以(suy)方程组的一般解为方程组的一般解为 31243529298921 1135270270TTXk第34页/共131页第三十五页，共132页。3、将矩阵、将矩阵(j zhn)对角化对角化 120230302A a=-1 2 0;-2 3 0;3 0 2; v,d=eig(a)v = 0 379/1257 379/1257 0 379/1257 379/1257 1 -379/419 -379/4

21、19 d =2 0 0 0 1 0 0 0 1 1VAVdA的特征值为的特征值为2，1，1 第35页/共131页第三十六页，共132页。4、用正交变换化二次型为标准、用正交变换化二次型为标准(biozhn)形形 2224123412131423243422 2222fxxxxx xx xx xx xx xx x 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

22、-0.6124 0 0.5000t = -0.0000 0 0 0 0 -0.0000 0 0 0 0 0 0 0 0 0 4.0000第36页/共131页第三十七页，共132页。 a=1 1 1 1；1 1 1 1；1 1 1 1；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 “*”表示(bios

23、h) 0 0 0 0 近似于零 0 0 0 4 FORMAT RAT Approximation by ratio of small integers.第37页/共131页第三十八页，共132页。4、用正交变换化二次型为标准、用正交变换化二次型为标准(biozhn)形形 2224123412131423243422 2222fxxxxx xx xx xx xx xx x结论结论(jiln)：作正交变换：作正交变换 112233440.08460.49280.70710.50000.08460.4928-0.70710.5000=-0.7815-0.373200.50000.6124-0.612

24、400.5000 xyxyxyxy则有则有 244fy第38页/共131页第三十九页，共132页。上机实验题一、基础型实验1、计算下列(xili)极限xeexxxsinlim0nnmmaxaxaxlimnxxx21lim111limxxe111limxxe第39页/共131页第四十页，共132页。2、计算下列(xili)导数(1)(2)(3)(4)1ln(2xxeeyxey1sin2212arcsinttyxxy 第40页/共131页第四十一页，共132页。一. 输入(shr)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.

25、 A-B; 3. A*B; 4. A*u; 5. 2A-3B; 6. A2+B2; 7. AB-BA。二. 求下列(xili)矩阵的逆阵并求其行列式的值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=2,1,-1;2,1,0;1,-1,1,B=1,-1,3;4,3,2;1,-2,5,XA=B; 3.A=

26、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.第41页/共131页第四十二页，共132页。 四. 将下列(xili)矩阵化为阶梯矩阵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. 五.求下列(xili)矩阵的秩 1.A=-5,6,-3;3,

27、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.第42页/共131页第四十三页，共132页。MATLAB典型函数含义MATLAB典型函数含义abs(x)求绝对值tan(x)正切值sqrt(x)求平方根值cot(x)余切值exp(x)指数运算atan(x)反正切值sin(x)正弦值acot(x)反余切值cos(x)余弦值log(x)自然对数asin(x)反正弦值Log10(x)常用对数acos(x)反余弦值 附录

28、：MATLAB软件(run jin)中部分常用函数表第43页/共131页第四十四页，共132页。 一、绘图功能第44页/共131页第四十五页，共132页。nx=0:pi/100:2*pi;ny=sin(x);nplot(x,y)第45页/共131页第四十六页，共132页。第46页/共131页第四十七页，共132页。【例3】 用不同线型和颜色重新绘制例2图形，其程序为：x=0:pi/100:2*pi;y1=sin(x);y2=cos(x);plot(x,y1,go,x,y2,b-.)其中参数(cnsh)go和b-.表示图形的颜色和线型。g表示绿色，o表示图形线型为圆圈；b表示蓝色，-.表示图形线

29、型为点划线。第47页/共131页第四十八页，共132页。第48页/共131页第四十九页，共132页。第49页/共131页第五十页，共132页。naxis （equal）两个坐标两个坐标因子设成相等因子设成相等naxis （off）关闭坐标系关闭坐标系统统naxis （on）显示坐标系显示坐标系统统第50页/共131页第五十一页，共132页。第51页/共131页第五十二页，共132页。【例5】 为正弦、余弦(yxin)曲线增加图例，其程序为：x=0:pi/100:2*pi;y1=sin(x);y2=cos(x);plot(x,y1,x,y2, -);legend(sin(x),cos(x);第5

30、2页/共131页第五十三页，共132页。阅读以下程序：阅读以下程序：x=-2:0.1:2; %产生产生(chnshng)横坐标横坐标x数组数组y=x.3-3*x; %计算由计算由y=x3-3x确定的纵坐标确定的纵坐标y数组数组plot(x,y) %绘图绘图grid on %给图形加上网格线给图形加上网格线axis equal %使使x,y轴单位刻度相等轴单位刻度相等第53页/共131页第五十四页，共132页。二、二、 subplot 并列并列(bngli)绘图绘图命令命令第54页/共131页第五十五页，共132页。第55页/共131页第五十六页，共132页。图形窗口的标题栏中，即图形图形窗口的

31、标题栏中，即图形窗口标题。用户可通过句柄激窗口标题。用户可通过句柄激活或关闭某图形窗口，而活或关闭某图形窗口，而axis、xlabel、title等许多命令也只对等许多命令也只对活动窗口有效。活动窗口有效。第56页/共131页第五十七页，共132页。1); nH2=figure; %创建第二个窗口并创建第二个窗口并返回句柄到变量返回句柄到变量H2nplot(x,z); %绘制绘制(huzh)图形并图形并设置有关属性设置有关属性ntitle(cos(x); axis (0 2*pi -1 1);nH3=figure; %同上同上nplot(x,t); title(tangent(x); axis

32、 (0 2*pi -40 40);nH4=figure; %同上同上nplot(x,ct); title(cotangent(x); axis (0 2*pi -40 40);第57页/共131页第五十八页，共132页。第58页/共131页第五十九页，共132页。线线naxis (0 2*pi -1 1); nlegend(cos,sin);nhold off %关闭图形关闭图形保持保持第59页/共131页第六十页，共132页。n a,b为为x的区间，的区间，c,d为为y的的区间。区间。n例：例：fplot(sin(x),0 2*pi,-+)nfplot(sin(x),cos(x),0 2*p

33、i,.) %同时绘制正弦、余弦曲线同时绘制正弦、余弦曲线第60页/共131页第六十一页，共132页。第61页/共131页第六十二页，共132页。2*pi)n%绘制隐函数绘制隐函数f(x,y)=0在在a,b与与c,d区间上的图形区间上的图形nezplot(4*x2+16*y2-3,-1 1 -1 1)n%绘制参数方程绘制参数方程x=sinx,y=cosx的图形的图形nezplot(sin(x),cos(x),0 2*pi)第62页/共131页第六十三页，共132页。第63页/共131页第六十四页，共132页。nx=0:0.01:2*piny=abs(1000*sin(4*x)+1nsemilog

34、y(x,y); %单对数单对数Y轴绘轴绘图命令图命令第64页/共131页第六十五页，共132页。ntitle(polar plot);第65页/共131页第六十六页，共132页。第66页/共131页第六十七页，共132页。第67页/共131页第六十八页，共132页。ny=0 0 1 1 0;nfill(x,y,y);%绘制并以黄色填充绘制并以黄色填充正方形图正方形图第68页/共131页第六十九页，共132页。出各种颜色。出各种颜色。第69页/共131页第七十页，共132页。第70页/共131页第七十一页，共132页。npolar 极坐标极坐标图图nbar 条形图条形图nloglog 双对数双对

35、数坐标图坐标图nsemilogx X轴轴为对数的坐标图为对数的坐标图nsemilogy Y轴轴为对数的坐标图为对数的坐标图nstairs 阶梯形阶梯形图图axis 设置设置(shzh)坐标轴坐标轴figure 创建图形窗口创建图形窗口grid 放置坐标网格线放置坐标网格线hold 保持当前图形窗口内容保持当前图形窗口内容subplot 创建子图创建子图title 放置图形标题放置图形标题xlabel 放置放置X轴坐标标记轴坐标标记ylabel 放置放置Y轴坐标标记轴坐标标记第71页/共131页第七十二页，共132页。阅读下面程序阅读下面程序(chngx):%绘制摆线绘制摆线:hold ont=

36、0:0.01:4*pi;for a=1:1:3 x=a*(t-sin(t); y=a*(1-cos(t); plot(x,y)end第72页/共131页第七十三页，共132页。h=3 2 1 0.5; %在曲线上取不同的点在曲线上取不同的点a=(exp(h)-1)./h;%计算连接点计算连接点M与与点与与点P的各条割线的各条割线(gxin)的斜率的斜率x=-1:0.1:3; %选定图形的自变量范围选定图形的自变量范围plot(x,exp(x),r);%作函数图形作函数图形hold on; %在图形上继续作图在图形上继续作图for i=1:4 plot(h(i),exp(h(i),w) %在图上

37、作出不同的点在图上作出不同的点 plot(x,a(i)*x+1) %作割线作割线(gxin)的图的图endaxis square %把所有图形放在一个正方形框内把所有图形放在一个正方形框内plot(x,x+1,g) %画出切线的图形画出切线的图形画出画出 在点在点P(0,1)处的切线及若干条割线处的切线及若干条割线(gxin),观察观察割线割线(gxin)的变化趋势的变化趋势,理解导数的定义及几何意义理解导数的定义及几何意义.xexf)(第73页/共131页第七十四页，共132页。n函数功能：以向量函数功能：以向量x，y，z为坐为坐标，绘制三维曲线。标，绘制三维曲线。第74页/共131页第七十

38、五页，共132页。第75页/共131页第七十六页，共132页。，组成了三维空间的网格点；组成了三维空间的网格点；c用用于控制网格点颜色。于控制网格点颜色。【例10】 下列程序绘制三维网格(wn )曲面图x=0:0.15:2*pi;y=0:0.15:2*pi;z=sin(y)*cos(x); %矩阵相乘mesh(x,y,z);第76页/共131页第七十七页，共132页。【例11】 下列程序绘制三维曲面(qmin)图形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)

39、,zlabel(z-label);title(3-D surf);第77页/共131页第七十八页，共132页。例绘制马鞍面的图形，并用平行截例绘制马鞍面的图形，并用平行截面法观察马鞍面的特点面法观察马鞍面的特点22zxyx=-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 onmesh(mx,my,mz)plot3(px,py,pz,r*)subplot(1,2,2)plot3(px,py,pz)第7

40、8页/共131页第七十九页，共132页。拟拟 合合2.2.拟合拟合(n h)(n h)的基的基本原理本原理1. 拟合(n h)问题引例第79页/共131页第八十页，共132页。拟拟 合合 问问 题题 引引 例例 1 1温度温度t(0C) 20.5 32.7 51.0 73.0 95.7电阻电阻R( ) 765 826 873 942 1032已知热敏电阻数据：已知热敏电阻数据：求求600C600C时的电阻时的电阻(dinz)R(dinz)R。 设 R=at+ba,b为待定系数(xsh)第80页/共131页第八十一页，共132页。拟拟 合合 问问 题题 引引 例例 2 2 t (h) 0.25

41、0.5 1 1.5 2 3 4 6 8c ( g/ml) 19.21 18.15 15.36 14.10 12.89 9.32 7.45 5.24 3.01已知一室模型快速静脉注射下的血药浓度数据已知一室模型快速静脉注射下的血药浓度数据(t=0注射注射300mg)求血药浓度随时间求血药浓度随时间(shjin)的变化规律的变化规律c(t).作半对数作半对数(du sh)坐标系坐标系(semilogy)下的图下的图形形为待定系数kcectckt,)(002468100101102MATLAB(aa1)第81页/共131页第八十二页，共132页。曲曲 线线 拟拟 合合 问问 题题 的的 提提 法法已

42、知一组（二维）数据，即平面上已知一组（二维）数据，即平面上 n n个点（个点（xi,yi) i=1,n, xi,yi) i=1,n, 寻求寻求一个函数（曲线）一个函数（曲线）y=f(x), y=f(x), 使使 f(x) f(x) 在某种准则下与所有数据点最为在某种准则下与所有数据点最为接近接近(jijn)(jijn)，即曲线拟合得最好。，即曲线拟合得最好。 +xyy=f(x)(xi,yi)i i 为点（为点（xi,yi) 与曲线与曲线(qxin) y=f(x) 的距离的距离第82页/共131页第八十三页，共132页。拟合拟合(n h)与插值的关系与插值的关系 函数(hnsh)插值与曲线拟合都

43、是要根据一组数据构造一个函数(hnsh)作为近似，由于近似的要求不同，二者的数学方法上是完全不同的。 实例：下面数据是某次实验所得，希望实例：下面数据是某次实验所得，希望(xwng)得到得到X和和 f之间的关系？之间的关系？x124791 21 31 51 7f1 .53 .96 .611 .71 5 .61 8 .81 9 .62 0 .62 1 .1MATLAB(cn)问题：问题：给定一批数据点，需确定满足特定要求的曲线或曲面解决方案：解决方案：若不要求曲线（面）通过所有数据点，而是要求它反映对象整体的变化趋势，这就是数据拟合数据拟合，又称曲线拟合或曲面拟合。若要求所求曲线（面）通过所给所

44、有数据点，就是插值问题插值问题；第83页/共131页第八十四页，共132页。最临近最临近(ln jn)插值、线性插值、样条插值与曲线拟合结果：插值、线性插值、样条插值与曲线拟合结果：第84页/共131页第八十五页，共132页。曲线拟合问题曲线拟合问题(wnt)最常用的解法最常用的解法线性最小二乘法的基本思路线性最小二乘法的基本思路第一步:先选定(xun dn)一组函数 r1(x), r2(x), rm(x), m0)k(0)模型模型(mxng)假设假设1. 1. 机体看作一个房室，室内血药浓度均匀机体看作一个房室，室内血药浓度均匀一室模型一室模型模型建立模型建立d/c(0) 3得：由假设-kc

45、dtdc 2得：由假设ktevdtc)( 在此，在此，d=300mg，t及及c（t）在某些点处的值见前表，需经）在某些点处的值见前表，需经拟合求出参数拟合求出参数k、v第108页/共131页第一百零九页，共132页。用线性最小二乘拟合用线性最小二乘拟合(n (n h)c(t)h)c(t)ktevdtc)()/ln(,ln21vdakacyktvdc)/ln(ln2/,121aedvakatayMATLAB(lihe1)计算结果：计算结果：)(02.15),/1 (2347. 0lvhkd=300;t=0.25 0.5 1 1.5 2 3 4 6 8;c=19.21 18.15 15.36 14

46、.10 12.89 9.32 7.45 5.24 3.01;y=log(c);a=polyfit(t,y,1)k=-a(1)v=d/exp(a(2)程序：程序：用非线性最小二用非线性最小二乘拟合乘拟合(n (n h)c(t)h)c(t)第109页/共131页第一百一十页，共132页。给药方案给药方案(fng n) 设计设计cc2c10t 设每次注射剂量D, 间隔时间 血药浓度c(t) 应c1 c(t) c2 初次(ch c)剂量D0 应加大,0DD给药方案记为：给药方案记为：kecc2112ln1cck2、)( ,1220ccDcD1、计算结果：计算结果：9 . 3, 3 .225, 5 .3

47、750DD)(4),(225),(3750hmgDmgD给药方案：给药方案：c1=10,c2=25k=0.2347v=15.02第110页/共131页第一百一十一页，共132页。故可制定故可制定(zhdng)给药方案：给药方案：)(4),(225),(3750hmgDmgD即即: 首次首次(shu c)注射注射375mg， 其余每次注射其余每次注射225mg， 注射的间隔时间为注射的间隔时间为4小时。小时。第111页/共131页第一百一十二页，共132页。估计水塔估计水塔(shut)的流量的流量2、解题、解题(ji t)思路思路3、算法、算法(sun f)设计设计与编程与编程1、问题、问题第1

48、12页/共131页第一百一十三页，共132页。 某居民区有一供居民用水的园柱形水塔，一般可以通过测量其水位来估计水的流量，但面临的困难是，当水塔水位下降到设定的最低水位时，水泵自动启动向水塔供水，到设定的最高水位时停止供水，这段时间无法(wf)测量水塔的水位和水泵的供水量通常水泵每天供水一两次，每次约两小时.水塔是一个高12.2米，直径17.4米的正园柱按照设计，水塔水位降至约8.2米时，水泵自动启动，水位升到约10.8米时水泵停止工作表1 是某一天的水位测量记录，试估计任何时刻（包括水泵正供水时）从水塔流出的水流量，及一天的总用水量第113页/共131页第一百一十四页，共132页。 表 1

49、水位测量记录 （符号/表示水泵启动）时刻(h)水位(cm)0 0.92 1.84 2.95 3.87 4.98 5.90 7.01 7.93 8.97968 948 931 913 898 881 869 852 839 822时刻(h)水位(cm)9.98 10.92 10.95 12.03 12.95 13.88 14.98 15.90 16.83 17.93/ / 1082 1050 1021 994 965 941 918 892时刻(h)水位(cm)19.04 19.96 20.84 22.01 22.96 23.88 24.99 25.91866 843 822 / / 1059

50、1035 1018第114页/共131页第一百一十五页，共132页。流量估计的解题流量估计的解题(ji t)思路思路拟合(n h)水位时间函数确定流量(liling)时间函数估计一天总用水量第115页/共131页第一百一十六页，共132页。 拟合水位拟合水位时间函数时间函数 测量记录看，一天有两个供水时段（以下称第测量记录看，一天有两个供水时段（以下称第1供供水时段和第水时段和第2供水时段），和供水时段），和3个水泵不工作时段（以个水泵不工作时段（以下称第下称第1时段时段t=0到到t=8.97，第，第2次时段次时段t=10.95到到t=20.84和第和第3时段时段t=23以后）对第以后）对第1

51、、2时段的测量数据直时段的测量数据直接分别作多项式拟合，得到水位函数为使拟合曲线接分别作多项式拟合，得到水位函数为使拟合曲线比较光滑，多项式次数不要太高，一般比较光滑，多项式次数不要太高，一般(ybn)在在36由于第由于第3时段只有时段只有3个测量记录，无法对这一时个测量记录，无法对这一时段的水位作出较好的拟合段的水位作出较好的拟合第116页/共131页第一百一十七页，共132页。 2、确定流量时间函数 对于(duy)第1、2时段只需将水位函数求导数即可，对于(duy)两个供水时段的流量，则用供水时段前后（水泵不工作时段）的流量拟合得到，并且将拟合得到的第2供水时段流量外推，将第3时段流量包含

52、在第2供水时段内第117页/共131页第一百一十八页，共132页。3、一天总用水量的估计 总用水量等于两个水泵不工作时段和两个供水时段用水量之和，它们(t men)都可以由流量对时间的积分得到。第118页/共131页第一百一十九页，共132页。算法算法(sun f)设计与设计与编程编程1、拟合、拟合(n h)第第1、2时段的水位，并导出流量时段的水位，并导出流量2、拟合、拟合(n h)供水时段的流量供水时段的流量3、估计一天总用水量估计一天总用水量4、流量及总用水量的检验、流量及总用水量的检验第119页/共131页第一百二十页，共132页。 1、拟合第1时段的水位，并导出流量 设t，h为已输入

53、的时刻和水位测量(cling)记录（水泵启动的4个时刻不输入），第1时段各时刻的流量可如下得：1） c1=polyfit（t（1：10），h（1：10），3）； %用3次多项式拟合第1时段水位，c1输出3次多项式的系数2）a1=polyder（c1）； % a1输出多项式（系数为c1）导数的系数 3）tp1=0：0.1：9； x1=-polyval（a1，tp1）； % x1输出多项式（系数为a1）在tp1点的函数值（取负后边为正值），即tp1时刻的流量 MATLAB(llgj1)4）流量(liling)函数为：1079.227173. 22356. 0)(2tttf第120页/共131页第一

54、百二十一页，共132页。 2、拟合第2时段的水位，并导出流量 设t，h为已输入(shr)的时刻和水位测量记录（水泵启动的4个时刻不输入(shr)），第2时段各时刻的流量可如下得：1） c2=polyfit(t(10.9:21),h(10.9:21),3)； %用3次多项式拟合第2时段水位，c2输出3次多项式的系数2） a2=polyder(c2)； % a2输出多项式（系数为c2）导数的系数 3）tp2=10.9:0.1:21; x2=-polyval(a2,tp2); % x2输出多项式（系数为a2）在tp2点的函数值（取负后边为正值），即tp2时刻的流量MATLAB(llgj2)4）流量(

55、liling)函数为：8313. 17512. 87529. 00186. 0)(23ttttf第121页/共131页第一百二十二页，共132页。 3、拟合供水时段的流量 在第1供水时段（t=911）之前(zhqin)（即第1时段）和之后（即第2时段）各取几点，其流量已经得到，用它们拟合第1供水时段的流量为使流量函数在t=9和t=11连续，我们简单地只取4个点，拟合3次多项式（即曲线必过这4个点），实现如下： xx1=-polyval（a1，8 9）； %取第1时段在t=8,9的流量 xx2=-polyval（a2，11 12）； %取第2时段在t=11,12的流量 xx12=xx1 xx2；

56、 c12=polyfit（8 9 11 12，xx12，3）； %拟合3次多项式 tp12=9：0.1：11； x12=polyval（c12，tp12）； % x12输出第1供水时段 各时刻的流量MATLAB(llgj3)拟合的流量拟合的流量(liling)函数为：函数为：078.3555879.737207. 3)(2tttf第122页/共131页第一百二十三页，共132页。 在第2供水时段之前取t=20，20.8两点的流水量，在该时刻之后（第3时段）仅有3个水位记录，我们用差分得到流量，然后用这4个数值拟合第2供水时段的流量如下： dt3=diff（t(22：24)）； %最后(zuhu

57、)3个时刻的两两之差 dh3=diff（h(22：24)）； %最后(zuhu)3个水位的两两之差 dht3=-dh3./dt3； %t(22)和t(23)的流量 t3=20 20.8 t(22) t(23)； xx3=-polyval(a2，t3(1：2)，dht3)； %取t3各时刻的流量 c3=polyfit（t3，xx3，3）； %拟合3次多项式 t3=20.8：0.1：24； x3=polyval（c3，tp3）；% x3输出第2供水时段 （外推至t=24）各时刻的流量MATLAB(llgj4)拟合的流量拟合的流量(liling)函数为：函数为：8283.913077. 71405. 0)(2tttf第123页/共131页第一百二十四页，共132页。 3、一天总用水量的估计 第1、2时段和第1、2供水时段流量的积分之和，就是一天总用水量虽然诸时段的流量已表为多项式函数，积分可以解析地算出，这里(zhl)仍用数值积分计算如下： y1=0.1*trapz(x1)； %第1时段用水量（仍按高 度计），0.1为积分步长 y2=0.1*trap