Matlab编程-第五章拟合与插值

第五章拟合与插值t(h)0.250.511.523468u(g/ml)19.2118.1515.3614.1012.899.327.455.243.01例5.1试验得快速静脉注射下的一组血药浓度数据求血药浓度随时间的变化规律u(t).(t=0时注射300mg，药物排除速率与血药浓度成正比)解：记血液容积v,t=0注射剂量d,血药浓度即为d/v.kudtduvdu)0(,kdtudu,dtkudu1lnCktuktCeuvdu)0(ktevdu如何确定常数k和v？一、用最小二乘法确定函数中的参数问题提出：对函数(曲线)y=f(x,k),k为待定参数，已测得一组数据，即曲线上的n个点(xi,yi)i=1,…,n,试确定参数k。(k=[k1,k2,…,km],mn)最小二乘准则：使n个点(xi,yi)与曲线y=f(x,k)的距离i的平方和最小：2211minmin[(,)]nniiiiifxkyMatlab中最小二乘拟合函数：(1)lsqcurvefit[a,J]=lsqcurvefit(fun,a0,x,y,lb,ub,options)fun:自定义的要拟合成的函数a0：参数的初值x，y：数据向量，lb，ub：a的上下限a：所求的参数J：J=min(sum(yi-f(a,xi))^2)最小平方误差(2)lsqnonlina=lsqnonlin(fun,a0,x,y)fun:自定义的要拟合成的函数a0：参数的初值x，y：数据向量，lb，ub：a的上下限a：所求的参数lsqnonlin可用以求向量值函数12(,)[(,),(,),...,(,)]Ttfxkfxkfxkfxk中的参数。例5.2薄膜渗透率问题某医用薄膜允许一种物质分子穿透，从高浓度溶液向低浓度溶液扩散，试制薄膜需测定穿透率。测定方法：用面积为S的薄膜将容器分为A、B两部分，体积为Va与Vb，分别注满不同浓度的溶液。该物质通过单位面积薄膜的速度与两侧溶液的浓度差成正比，比例系数为K称为穿透率。定时测量一侧的浓度，以确定K。已知Va=Vb=1000立方厘米，S=10平方厘米，以及B侧浓度表：时间t单位：s，浓度Cb单位：10-5mg/cm3t1002003004005006007008009001000Cb454499545576596615625638645650假设：（1）两侧各自浓度均匀（2）分子总是从高浓度溶液向低浓度溶液扩散符号：Ca(t),Cb(t)表示t时刻两侧的浓度。分析：考虑时段：],[tttA侧物质质量的变化：)()(tCVttCVaaaattCtCSKab)]()([ttCttCaa)()()]()([tCtCVSKabattCttCaat)()(lim0)]()([tCtCVSKaba)]()([)(tCtCVSKtCabaa整个容器物质总质量不变：)0()0()()(bbaabbaaCVCVtCVtCV由Va=Vb：)()0()0()(tCCCtCbbaa.......(*)(*)式化为：)]0()0()(2[)(bababCCtCVSKtC)]0()0([)(2)(baababCCVSKtCVSKtC)]0()0([)(2)(baababCCVSKtCVSKtCtVSKbaabatVSKtVSKbaaaeCCVSKtCVSKeetC222)]0()0([)(2)(tVSKbaatVSKbaaeCCVSKetC22)]0()0([)(CeCCetCtVSKbatVSKbaa22)]0()0([21)(CeCCetCtVSKbatVSKbaa22)]0()0([21)(CCCCbab)]0()0([21)0()]0()0([21abCCCtVSKabbabaeCCCCtC2)]0()0([21)]0()0([21)(问题化为：由实验数据来识别参数tVSKabbabaeCCCCtC2)]0()0([21)]0()0([21)(记：)]0()0([21)],0()0([21abbaCCbCCaKtbbeatC02.0)(t1002003004005006007008009001000Cb454499545576596615625638645650KtbbeatC02.0)(functionf=fun502(x,t)f=x(1)+x(2)*exp(-0.02*x(3)*t);t=100:100:1000;cb=10^(-5)*[454,499,545,576,596,615,625,638,645,650];x0=[0.05,-0.005,0.05];x=lsqcurvefit('fun502',x0,t,cb)t1=100:20:1000;y=fun502(x,t1);plot(t,cb,'r*',t1,y)KtbbeatC02.0)(得：a=0.0066,b=-0.0030k=0.1486100200300400500600700800900100044.555.566.57x10-3例5.3在某次阻尼振荡实验中测得18组数据点，试确定其振动方程。t：00.41.222.83.64.45.26y:10.850.29-0.27-0.53-0.4-0.120.170.28t：7.289.210.411.612.413.614.415y:0.15-0.03-0.15-0.070.060.080.032-0.015-0.02由物理背景知：wtektay)cos(clc,clearx=[0,0.4,1.2,2,2.8,3.6,4.4,5.2,6,7.2,8,9.2,10.4,11.6,12.4,13.6,14.4,15]';y=[1,0.85,0.29,-0.27,-0.53,-0.4,-0.12,0.17,0.28,0.15,-0.03,-0.15,-0.07,0.059,0.08,0.032,-0.015,-0.02]';f=fittype('a*cos(k*x)*exp(w*x)','independent','x','coeff',{'a','k','w'});cfun=fit(x,y,f)xi=0:0.1:20;yi=cfun(xi);plot(x,y,'r*',xi’,yi’,'b-')由物理背景知：wtektay)cos(运行结果：Generalmodel:cfun(x)=a*cos(k*x)*exp(w*x)Coefficients(with95%confidencebounds):a=0.9988(0.9837,1.014)k=1.001(0.9959,1.006)w=-0.2067(-0.2132,-0.2003)tety2067.0-)001.1cos(9988.0051015-0.6-0.4-0.200.20.40.60.81二、曲线拟合：曲线拟合问题的提法：已知曲线上n个点(xi,yi),(i=1,2,…,n),xi互不相同，寻求一个函数y=F(x),使y=F(x)在某种准则下与所有数据点最为接近，即曲线拟合得最好。最小二乘法是解决曲线拟合最常用的方法，基本思路是:mfff,,,10设是m+1个线性无关的连续函数，称为基本拟合函数，取mmfcfcfcxF1100)(mccc,,,10为待定系数，mn，niiiyxF12])([min拟合准则是:按此方法得到的函数mmfcfcfcxF1100)(称为数据集(xi,yi),i=1,2,…,n的最小二乘拟合函数。衡量拟合情况优劣的量之一是：222)-())(-(1---1yyxFynmmRiii其中miiymy011R2的值可作为拟合质量的指标。若它接近1，表示F(x)能较好的拟合该数据集。Matlab中的多项式拟合函数：polyfit格式：a=polyfit(x,y,n)功能：对一组点(x,y)进行n次多项式拟合。x,y为要拟合的数据，n为拟合多项式的次数，a为拟合多项式系数，},,,,1{},,,{210mmxxxfff若取：即用m次多项式拟合给定数据，称多项式拟合。0111axaxaxannnn],,,,[011aaaaann拟合多项式在x处的值可用polyval(a,x)计算。例5.4某通信公司在一次施工中，需要在一段水面宽20米的河底沿直线走向铺设一条河底光缆，A点表示光缆入水处，B点表示光缆出水处，C、D点表示河底的光缆的两个触底点，在铺设前对河底地形做了测量，测得一组等分点位置xi处的水深hi，x0位于点A处，x20位于点B处，数据如表（单位：米）试估算通过这条河沟所需光缆长度的近似值。ABCDx0x1x2x3x4x5x6x7x8x9x109.28.918.027.958.119.0510.1211.0812.2513.2613.45x11x12x13x14x15x16x17x18x19x2012.6211.2810.239.257.918.058.869.8510.8110.96河底深度测量数据（单位：米）得出对河底光缆长度的估计s，所需光缆为s+y(0)+y(20)对数据求出合适次数的拟合多项式p(x)，由弧长公式：dxxps2002))((1x=0:1:20;y=[-9.2,-8.91,-8.02,-7.95,-8.11,-9.05,-10.12,-11.08,-12.25,-13.26,-13.45,-12.62,-11.28,-10.23,-9.25,-7.91,-8.05,-8.86,-9.85,-10.81,-10.96];polytool(x,y,3)n=input('拟合次数n=')formatlongp=polyfit(x,y,n)dp=polyder(p)x1=0:0.01:20;y1=sqrt(1+polyval(dp,x1).^2);s=trapz(x1,y1)-y(1)-y(end)6次多项式拟合较为合适。光缆总长：46.38m三、一维插值已知函数y=f(x)在区间[a,b]上的n+1个不同点的函数值为若存在一个简单函数F(x),使称F(x)为f(x)在区间[a,b]上的插值函数，称(xi,yi)为插值节点。nxxx10nyyy,,,10niyxFii,,2,1,0,)(当函数的表示式未知或非常复杂时，可通过插值的方法研究函数。若F(x)为多项式，称为多项式插值(或代数插值)；常用的代数插值方法有：拉格朗日插值，牛顿插值等。n次代数插值：已知函数f(x)在n+1个点x0,x1,…,xn处的函数值为y0,y1,…,yn，求一n次多项式函数Pn(x)，使其满足：Pn(xi)=yi,(i=0,1,…,n).若Pn(x)按下述方式构造，称为拉格朗日插值，称为拉格朗日插值基函数．其中Li(x)为n次多项式：0()()nniiiPxLxy01110111()()()()()()()()()()()iiniiiiiiiinxxxxxxxxxxLxxxxxxxxxxx特别地:（1）已知两个节点时，得线性插值多项式:（2）已知三个节点时，得抛物插值多项式:01100110xxxxLxyyxxxx020112201010210122021xxxxxxxxxxxxLxyyyxxxxxxxxxxxx已知n+1个节点时，可得n次拉格朗日插值多项式。在区间[-1,1]上的代数插值22511)(xxf当插值次数过高时，常常会出现龙格现象，因此在数据建模中不建议使用过高次的代数插值。若F(x)为分段函数，称为分段插值；插值中使用较多的是分段线性插值和三次样条插值。1.分段线性插值设函数y=