6医用薄膜渗透率的确定―数据拟合

下一页上一页主页实验名称：曲线拟合方法建模与计算实验目的：熟悉并掌握数据拟合的基本原理、方法及MATLAB实现实验任务：1、医用薄膜渗透率的确定——建模与计算2、教材第87页操练一实验五下一页上一页主页实验目的[1]了解最小二乘拟合的原理，掌握用MATLAB作线性最小二乘拟合的方法。[2]通过实例学习如何用拟合方法解决实际问题，注意与插值方法的区别；通过实例理解参数辨识问题的几种方法。[3]鼓励不囿于固定的模式或秩序，灵活调整思路，突破思维的呆板性，找到打破常规的解法。并在文献检索、动手和动脑等方面得到锻炼，树立创新意识。实验内容范例：薄膜渗透率的测定医用薄膜的渗透率拟合的基本原理拟合函数的选取用MATLAB作拟合计算下一页上一页主页某种医用薄膜在试制时需测定其被物质分子穿透的能力。测定方法：用面积为S的薄膜将容器分成两部分，在两部分中分别注满该物质的两种不同浓度的溶液。此时该物质分子就会从一侧向另一侧扩散。平均每单位时间通过单位面积薄膜的物质分子量与膜两侧溶液的浓度差成正比，比例系数K称为渗透率。定时测量容器中薄膜某一侧的溶液浓度，以此确定K。VAVBS问题背景返回医用薄膜的渗透率下一页上一页主页3）薄膜是双向同性的即物质从膜的任何一侧向另一侧渗透的性能是相同的。假设1）薄膜两侧的溶液始终是均匀的；2）在单位时间内通过单位面积薄膜的物质分子量与膜两侧溶液的浓度差成正比。返回下一页上一页主页第一步：通过机理分析确定容器一侧的浓度随时间的变化规律：C=CB(t)第二步：利用实验数据(tj,CB(tj))，求出其中的未知参数，包括渗透率K.解决问题的思路下一页上一页主页考察时段[t，t+Δt]薄膜的一侧容器中该物质质量的变化。1）以容器A侧为例，在时段[t，t+Δt]物质质量的增量为：求浓度随时间的变化规律由假设2,在时段[t，t+Δt]，从B侧渗透至A侧的该物质的质量为：下一页上一页主页于是有：两边除以Δt，并令Δt→0取极限再稍加整理即得：（1）求浓度随时间的变化规律下一页上一页主页分别表示在初始时刻两侧溶其中2)注意到任意时刻，整个溶液中含有该物质的质量,与初始时刻该物质的含量相同，因此液的浓度.求浓度随时间的变化规律下一页上一页主页加上初值条件：从而：返回代入式（1）：解之得：求浓度随时间的变化规律下一页上一页主页问题归结为利用CB在时刻tj的测量数据Cj(j=1,2,...,N)来确定（辨识）K和.待定参数的确定Tj(秒)1002003004005006007008009001000Cj(×103)4.544.995.355.655.906.106.266.396.506.59容器的B部分溶液浓度的测量数据下一页上一页主页曲线拟合的基本原理已知一组（二维）数据，即平面上n个点（xi,yi)i=1,…n,寻求一个函数（曲线）y=f(x),使f(x)在某种准则下与所有数据点最为接近，即曲线拟合得最好。+++++++++xyy=f(x)(xi,yi)ii为点（xi,yi)与曲线y=f(x)的距离下一页上一页主页问题的数学模型步骤：1)选定一类函数f(x，a1，a2，…，am)（1）其中a1,a2,…am为待定常数。2)确定参数a1,a2,…am，准则(最小二乘准则）：使n个点（xi,yi)与曲线y=f(x，a1，a2，…，am)的距离i的平方和最小。下一页上一页主页)2(]),,,([),,(2111221iniminiimyaaxfaaaJ记问题归结为:求a1,a2,…am使J(a1,a2,…am)最小。这样的拟合称为最小二乘拟合。问题的数学模型下一页上一页主页特别，若选定一组函数r1(x),r2(x),…rm(x),mn,令f(x,a1,a2,…am)=a1r1(x)+a2r2(x)+…+amrm(x)其中a1,a2,…am为待定系数。问题的解法——线性最小二乘法当f具有上述形式，且按最小二乘准则确定待定系数时，称这样的拟合为线性最小二乘拟合，称对应的函数f(x,a1,a2,…am)为已知数据点的线性最小二乘拟合函数。下一页上一页主页最小二乘拟合函数的选取+++++++++++++++f=a1+a2x1.将数据(xi,yi)i=1,…n作图，通过直观判断确定f：f=a1+a2x+a3x2f=a1+a2x+a3x2下一页上一页主页2.通过机理分析建立数学模型来确定f。+++++++++++++++f=a1+a2/xf=a1exp(a2x)f=a1exp(a2x)最小二乘拟合函数的选取返回下一页上一页主页1.多项式拟合：作多项式f(x)=a1xm+…+amx+am+1函数拟合,可利用已有程序polyfit,其调用格式为:a=polyfit(x,y,m)用MATLAB作最小二乘拟合数据点拟合多项式次数系数下一页上一页主页例.已知热敏电阻数据：温度t(0C)20.532.751.073.095.7电阻R()7658268739421032求电阻R随温度t的变化规律。2040608010070080090010001100用MATLAB作多项式拟合下一页上一页主页2.用命令polyfit(x,y,m)作最小二乘拟合3.编写MATLAB程序dianzu1.m，并运行得到：a1=3.3940,a2=702.49181.选取拟合函数R=a1t+a2用MATLAB作多项式拟合下一页上一页主页t=[20.532.5517395.7];r=[7658268739421032];aa=polyfit(t,r,1);a=aa(1)b=aa(2)y=polyval(aa,t);plot(t,r,'k+',t,y,'r')用MATLAB作多项式拟合下一页上一页主页2.曲线拟合：(1)作一般的最小二乘曲线拟合，可利用已有程序lsqcurvefit,其调用格式为：c＝lsqcurvefit(f,c0,x,y)用MATLAB作曲线拟合数据点c的初值函数M文件拟合函数y=f(a,x)的函数M—文件下一页上一页主页xdata=[0:.1:2]ydata=[5.89553.56392.51731.97901.89901.39381.13591.00961.03430.84350.68560.61000.53920.39460.39030.54740.34590.13700.22110.17040.2636]例：用函数y(x)=a1*exp(-a2*x)+a3*exp(-a4*x)拟合下列数据点：1.用命令lsqcurvefit(fun,a0,x,y)用MATLAB作曲线拟合下一页上一页主页fun=inline('a(1)*exp(-a(2)*x)+a(3)*exp(-a(4)*x)','a','x');a0=[1,1,1,0];xdata=0:.1:2;ydata=[5.89553.56392.51731.97901.89901.39381.13591.00961.03430.84350.68560.61000.53920.39460.39030.54740.34590.13700.22110.17040.2636];a=lsqcurvefit(fun,a0,xdata,ydata)xi=linspace(0,2,200);yi=fun(a,xi);plot(xdata,ydata,'ro',xi,yi)xlabel('x'),ylabel('y=f(x)'),title('nonlinearcurvefitting')用MATLAB作曲线拟合返回下一页上一页主页调用格式：[c,r]=lsqcurvefit(fun,c0,x,y)除了输出拟合参数c外，还可以输出残差。(2)、基于Gauss-Newton算法的最小二乘拟合beta=nlinfit(x,y,’fun’,beta0)其中，fun为自定义的非线性拟合模型，beta0为拟合参数初始值，输出值beta为拟合模型中的参数。下一页上一页主页[beta，r]=nlinfit(x,y,’fun’,beta0)调用格式：输出的r为此参数下残差构成的向量例：《数学模型》P29的拟合clearformatlongfun=inline('p(1)*n.^2+p(2).*n','p','n');t=[0,20,40,60,80,100,120,140,160,184]';n=[0,1141,2019,2760,3413,4004,4545,5051,5525,6061]';p0=[0,0];p3=lsqcurvefit(fun,p0,n,t)[p4,r]=lsqcurvefit(fun,p0,n,t)%p3=nlinfit(n,t,fun,p0)%[p3,r]=nlinfit(n,t,fun,p0)下一页上一页主页用拟合方法确定医用薄膜渗透率前面得到的容器B侧的浓度函数为：问题归结为利用曲线拟合方法和测量数据(tj,Cj)(j=1,2,...,N)来确定（辨识）CB中的参数K和。下一页上一页主页引入从而用函数CB(t)来拟合所给的实验数据，从而估计出其中的参数a,b,K。若将其代入上式有：拟合函数化简：下一页上一页主页用MATLAB软件进行计算.1）编写函数M-文件nongdu.mfunctionf=nongdu(x,tdata)f=x(1)+x(2)*exp(-0.02*x(3)*tdata);其中x(1)=a;x(2)=b;x(3)=k;2)编写M文件(bomo.m)tdata=linspace(100,1000,10);cdata=[4.544.995.355.655.906.10...6.266.396.506.59];x0=[0.2,0.05,0.05];x=lsqcurvefit(‘nongdu’,x0,tdata,cdata)编写MATLAB程序:下一页上一页主页3)输出结果:x=0.007-0.0030.1012即k=0.1012,a=0.007,b=-0.003,结果:进一步求得：返回下一页上一页主页结果及误差分析几种方法得出的结果及相应的误差总结于下表，误差为计算数据与实验数据之差的平方和。方法abk误差E函数拟合6.9852-2.99410.10123.5604e-007％作出拟合函数曲线和数据点的图形tdata=linspace(100,1000,10);cdata=1e-05.*[454499535565590...　　　　610626639650659];x0=[0.2,0.05,0.05];[x,resnorm,residual]=lsqcurvefit('nongdu',x0,tdata,cdata)t=linspace(100,1000,100);c=nongdu(x,t);plot(tdata,cdata,'o',t,c)％求拟合在各节点的误差平方和c1=nongdu(x,tdata);e=c1-cdata;e1=sum(e.*e)下一页上一页主页函数拟合法的拟合效果返回下一页上一页主页clearformatlongnlinfun=inline('y=x(1)+x(2)*exp(-0.02*x(3)*t)','x','t');t=[100:100:1000]';c=10^(-3)*[4.54,4.99,5.35,5.65,5.90,6.10,6.26,6.39,6.50,6.59]';p0=[0.2,0.05,0.05];[p3,r]=nlinfit(t,c,'nlinfun',p0)%norm(feval('nlinfun',p3,t)-c)^2norm(r)下一页上一页主页例：静脉注射的给药方案问题背景一种新药用于临床之前，必须设计给药方案。在快速静脉注射下，所谓给药方案是指，每次注射量多大，间隔时间多长。药物进入肌体后随血液输送到全身，在这过程中不断被吸收、分解、代谢，最终排出体外。药物向体外排出的速率与血药浓度成正比。单位体积血液中的药物含量，称血药浓度。临床上，每种药物有一个最小有效浓度c1和最大治疗浓度c2。2014-9-2535