分子穿透能力的测定

分子穿透能力的测定摘要:通过对问题的分析，根据质量守恒，利用微分方程模型，得到了关于浓度低一侧浓度对时间的微分方程模型。通过求解参数和简化以后确定了浓度与时间的指数关系。最后运用Matlab的最小二乘法进行拟合，求得参数，从而得到渗透率K。最后通过验证和检验发现求得结果与实验数据相吻合。关键词：渗透率质量守恒微分方程模型最小二乘法拟合Matlab参数估计一，问题重述:某种医用薄膜有允许一种物质的分子穿透它从高浓度的溶液向低浓度的溶液扩散的功能，在测试时需测定薄膜被这种分子穿透的能力。测定方法如下：用面积10cm2的薄膜将分成体积分别为100cm3和100cm3的两部分，在两部分中分别注满该物质的两种不同浓度的溶液。此时该物质分子就会从高浓度溶液穿过薄膜向低浓度溶液中扩散。通过单位面积膜分子扩散的速度与膜两侧溶液的浓度差成正比，比例系数K表征了薄膜被该物质分子穿透的能力，称为渗透率。定时测量容器中薄膜某一侧的溶液浓度值，以此确定K的数值。对容器一侧溶液浓度的测试结果如下：tj1002003004005006007008009001000Cj(10-3mg/cm3)4.544.995.355.655.906.106.266.396.506.59试建立一个较好的数学模型并给出相应的算法和程序。二，模型假设：1)薄膜两侧的溶液始终是均匀的，即在任何时刻薄膜两侧的每一处溶液的溶度都是相同的。2)当两侧浓度不一致时，物质的分子穿透薄膜总是从高浓度溶液向低浓度溶液扩散。3)通过单位面积膜分子扩散的速度与膜两侧溶液的浓度差成正比。4)薄膜是双向性的，即物质从膜的任何一侧向另一侧渗透的性能是相同的。三，符号说明：1)CA(t),CB(t)表示t时刻薄膜两侧溶液的浓度。2)aA，aB表示初始时刻薄膜两侧溶液的浓度（单位：mg/cm3）。3)K表示渗透率。4)VA,VB表示由薄膜阻隔的容器两侧的体积。四，问题的分析：1.通过机理分析确定容器一侧的浓度随时间的变化规律。2.利用实验数据（tj，Cj(tj)），求出未知的参数。3.最终进行检验。五，模型的建立：考察时段[t,t+Δt]薄膜两侧容器中该物质质量的变化。以容器A为例，在该时段物质质量曾加为：CAVA(t+Δt)-CAVA(t)另一方面有渗透率的定义可知，从B侧渗透至A侧的物质质量为:SK（CB−CA）Δt由质量守恒定律知两者相等，于是有：CAVA(t+Δt)-CAVA(t)=SK（CB−CA）Δt两边同时除以Δt，令Δt→0取极限，整理得：dCAdt=SKVA(CB−CA)①注意到整个容器的溶液中含有该物质的质量应该不变，即有下式成立：CAVA(t)+CBVB(t)=VAaA+VBaB②CA(t)=aA+VBVAaB--VBVACB(t)③将③式代入①式得：dCBdt+SK(1VA-+1VB)CB=SK(aBVA-+aAVB)初始条件CB(0)=aB。六，问题求解：1.函数拟合法前面得到的模型是一个带初值的一阶线性微分方程，解之得：CB(t)=VAaA+VBaBVA+VB+VA（aA−aB）VA+VBe−SK(1VA−+1VB)t问题归结为利用CB在时刻tj的测量数据Cj（j=1，2，…，N）来辨识K和aA，aB。令a=VAaA+VBaBVA+VB，b=VA（aA−aB）VA+VB则，CB(t)=a+be−SK(1VA−+1VB)t将VA=VB=100cm3，S=10cm2代入上式有：CB(t)=a+be−0.2Kt用函数CB(t)来拟合所给的实验数据，从而估计出其中的参数a，b，K。（1）编写M文件（nongdu.m）functionf=nongdu(x,tdata)f=x(1)+x(2)*exp(-0.2*x(3)*tdata);%ÆäÖÐx(1)=a;x(2)=b;x(3)=k;（2）编写主程序tdata=linspace(100,1000,10);cdata=[4.544.995.355.655.906.106.266.396.506.59];x0=[0.2,0.05,0.05];x=lsqcurvefit('nongdu',x0,tdata,cdata)（3）输出结果x=5.82701.92296.7010即表示结果K=6.7010，a=5.8270，b=1.9229进而求得aA=7.7499mg/cm3；aB=2.9041mg/cm32.非线性规划法利用CB在时刻tj的测量数据Cj（j=1，2，…，N）来辨识K和aA，aB。问题可转化为求函数：E(K,aA，aB)=∑(CB(tj)−𝐶j)2Nj=1⇒min即求函数E(K,aA，aB)=∑（a+be−0.2Ktj−𝐶j）2的最小值点（K,a,b）。下面用Matlab软件进行计算。（1）编写M文件（curvefun.m）functionf=curvefun(x,tdata)f=x(1)+x(2)*exp(-0.2*x(3)*tdata);%x(1)=a;x(2)=b;x(3)=k;（2）编写主程序tdata=linspace(100,1000,10);cdata=1e-03.*[4.544.995.355.655.906.106.266.396.506.59];x0=[0.2,0.05,0.05];x=lsqcurvefit('curvefun',x0,tdata,cdata)f=curvefun(x,tdata)plot(tdata,cdata,'o',tdata,f,'r-')（3）输出结果x=5.82701.92296.7010f=5.82705.82705.82705.82705.82705.82705.82705.82705.82705.8270即表示结果K=6.7010，a=5.8270，b=1.9229进而求得aA=7.7499mg/cm3；aB=2.9041mg/cm33.导函数拟合法前面得到的微分方程化简为：dCBdt=K[0.1(aA+aB)−0.2CB]令h=aA+aB上式变为：dCBdt=K[0.1h−0.2CB]这可以看作dCBdt随CB的变化规律。若知道一组数据[Cj，(dCBdt)j]（j=1，2，…，N）,则可以用最小二乘拟合的方法来求出函数dCBdt中的未知参数K，h。即为求参数K，h使得下列误差函数最小：J(K,h)=∑[(dCBdt)j−K(0.1h−0.2𝐶j)]2Nj=1该问题等价于用函数f(K，h，CB)=K（0.1h−0.2CB）来拟合数据[Cj，(dCBdt)j]（j=1，2，…，N）（1）编写M文件（baomof.m）functionf=baomof(x,cdata)f=x(1)*(0.1*x(2)-0.2*cdata)%x(1)=K;x(2)=h（2）编写主程序%求数据点[Cj，(dCBdt)j]（j=1，2，…，N）tdata=linspace(100,1000,10);cdata=[4.544.995.355.655.906.106.266.396.506.59];[d,ifail]=e01bef(tdata,cdata);[cj,dcj]=e01bgf(tdata,cdata,d,tdata);%作函数拟合x0=[0.2,0.1];x=lsqcurvefit('baomof',x0,cdata,dcj')（3）输出结果x=6.701010.6540即表示结果K=6.7010，h=10.65404.线性化迭代法前面带初始条件的一阶线性微分方程的解为：CB(t)=a+be−0.2Kt其中a=0.5(aA+aB)，b=(aA−aB)如果得到了参数K的一个较好的近似值K∗，则将CB(t)关于K在K∗处展开，略去ΔK的二次及以上的项，得CB(t)的一个近似式：CB(t)≈a+be−0.2K∗t+0.2b·ΔK·e−0.2K∗t·t通过极小化E(a，b，ΔK)=∑[𝐶j−(a+be−0.2Ktj+de−0.2K∗tj)]2Nj=1；其中d=0.2b·ΔK。确定a，b，d，再由ΔK=d/0.2b得到K∗的修正值ΔK。K∗←K∗·ΔK，得到K的一个新的近似值，用同样的方法再求新的修正值ΔK，这个过程可以不断重复，知道修正值足够小为止。