分子穿透能力的测定方案设计

数 学 建 模

第二次培训

分子穿透能力的测定方案设计

摘要

本文研究薄膜被物质分子穿透的能力，即求渗透率K。通过对问题的分析，根据质量守恒定律，得到了浓度低的一侧溶液的浓度与时间的微分方程模型。建模忽略了实验过程中溶液体积的改变和溶液溶质分布不均等的影响。运用1stOpt软件和Matlab软件编程拟合，求得渗透率K。最后对拟合结果进行检验。

关键词：渗透率 微分方程模型 质量守恒 参数估计

一、问题重述

某种医用薄膜有允许一种物质的分子穿透它从高浓度的溶液向低浓度的溶液扩散的功能，在测试时需测定薄膜被这种分子穿透的能力。测定方法如下：用面积10cm2的薄膜将分成体积分别为100cm3和100cm3的两部分，在两部分中分别注满该物质的两种不同浓度的溶液。此时该物质分子就会从高浓度溶液穿过薄膜向低浓度溶液中扩散。通过单位面积膜分子扩散的速度与膜两侧溶液的浓度差成正比，比例系数K表征了薄膜被该物质分子穿透的能力，称为渗透率。定时测量容器中薄膜某一侧的溶液浓度值，以此确定K的数值。

对容器一侧溶液浓度的测试结果如下：

tj 100 200 300 400 500 600 700 800 900 1000 Cj(10-3mg/cm3) 4.54 4.99 5.35 5.65 5.90 6.10 6.26 6.39 6.50 6.59 问题: 建立一个数学模型测定薄膜的渗透率并给出相应的算法和程序。

二、问题分析

渗透率K表征了薄膜被该物质分子穿透的能力，是薄膜的固定参数。该物质分子从高浓度溶液穿过薄膜向低浓度溶液中扩散，扩散速度与薄膜两侧溶液的浓度差成正比，比例系数为K。

在实验过程中，两部分溶液中溶质的总质量是不变的。设薄膜一侧的溶液为A，另一侧溶液为B；假设A溶液的浓度大于B溶液的浓度，则A溶液溶质质量的减少量等于B溶液溶质质量的增加量。

三、模型假设

1、薄膜是均匀的，即各单位面积薄膜分子扩散速度是相同的。 2、溶质在溶液中分布是均匀的。

3、在测试过程中忽略两部分溶液体积的改变，溶液体积始终等于初始体积100 cm3 。

4、两溶液中水分子的流动是双向的。

四、符号说明

序号 1 2 3 4 5 6 8 9 10

符号定义

符号说明 渗透率 A溶液的体积 B溶液的体积 A溶液的初始浓度 B溶液的初始浓度 实验进行的时间

K VA VB CA0 CB0 t CA(t) CB(t) S t时刻A溶液的浓度 t时刻B溶液的浓度

薄膜的面积

五、模型建立

在[t，t+]的时间内，溶液B中溶质质量的增加量为：

CB(tt)VBCBtVB在[t，t+ ①

]的时间内，穿过薄膜的溶质的质量为：

KSCAtCBtt ②

分析易知① = ②，即

CB(tt)VBCBtVB KSCAtCBt t ③

③等式两边同除以

得：

CB(tt)CBtKSCAtCBt ④ tVB当

0时，④可化为：

dCB(t)KSCAtCBt ⑤ dtVB在整个实验过程中两溶液中溶质的总质量是不变的，即

CA0VA CB0VB CAtVA CBtVB

化简得：

CB0CBtCAtCA0VB ⑥

VA将⑥代入⑤得

CtCBtCCdCB(t)KS(B)KS(A0B0) ⑦ dtVBVAVBVA等式⑦是一个一阶线性微分方程，其初始条件为：CB(0)=CB0 解得：

CA0VACB0VBCB0VACA0VA-SK(VAVB)t CBteVAVBVAVB11令 a =

CA0VACB0VBCVCA0VA b =B0A

VAVBVAVB-SK(11)tVAVB则 CBtabe ⑧

将VA = VB =100cm3 S=10cm2 带入⑧得：

CBtabe-0.2Kt ⑨

六、模型求解

由⑨可知CB(t)和t的关系，用该关系式拟合题目中给的定时测定的容器中薄膜一侧的溶液浓度值的实验数据，由此可确定a ,b , K的值。 方法一：

运用1stOp软件进行非线性拟合解得：

a = 0.00698504037105131 b = -0.00299407530539479 K = 0.0101227382188305

容器一侧的溶液浓度离散点与1stOpt拟合结果图见图6.1。

图6.1容器一侧的溶液浓度离散点与1stOpt拟合结果图

方法二:

运用Matlab进行非线性拟合解得： a = 0.006998(0.006978, 0.007019) b = -0.002999(-0.003014, -0.002985) K = 0.01(0.009821, 0.01018)

容器一侧的溶液浓度离散点与Matlab拟合结果图见图6.2。

图6.2 容器一侧的溶液浓度离散点与Matlab拟合结果图

七、结果分析

运用1stOpt软件对拟合的结果进行检验，检验的结果如下：

\"分子穿透能力测试\" 迭代数: 25

计算用时(时:分:秒:微秒): 00:00:00:563

优化算法: 麦夸特法(Levenberg-Marquardt) + 通用全局优化法 计算结束原因: 达到收敛判断标准 均方差(RMSE): 2.37761568253082E-6 残差平方和(SSE): 5.65305633381648E-11 相关系数(R): 0.999993326098118

相关系数之平方(R^2): 0.999986652240777 决定系数(DC): 0.999986652240777

卡方系数(Chi-Square): 4.94593781084979E-9 F统计(F-Statistic): 599343.536570763

参数 最佳估算 ----------

-------------

a 0.00698504037105131 b -0.00299407530539479

k 0.0101227382188305

====== 结果输出 =====

No 实测值y 计算值y 1 2 3 4 5 6 7 8 9

0.00454 0.0045397 0.00499 0.0049879 0.00535 0.0053539 0.00565 0.0056529 0.0059 0.0058970 0.0061 0.0060964 0.00626 0.0062593 0.00639 0.0063923 0.0065 0.0065009

10 0.00659 0.0065897

其中，均方差（RMSE）、残差平方和（SSE）、卡方系数（Chi-Square）的值越接近于0表示拟合的曲线和真实数据越接近；

相关系数（R）、相关系数平方（R2）、决定系数（DC）的值越接近1表示拟合程度越高。

运用Matlab软件对拟合的结果进行检验，检验的结果如下：

f = General model: f(t) = a+b*exp(-0.2*K*t)

Coefficients (with 95% confidence bounds): a = 0.006998(0.006978, 0.007019) b = -0.002999(-0.003014, -0.002985) K = 0.01(0.009821, 0.01018) goodness =struct with fields: sse: 9.0320e-11 rsquare: 1.0000 dfe: 7

adjrsquare: 1.0000 rmse: 3.5921e-06

其中，Sse表示误差平方和,值越小，表示拟合的曲线和真实数据越接近；

Rsquare表示确定系数,值越接近于1越好； dfe表示自由度即能够自由取值的变量个数 ；

Adjrsquare表示调整后的确定系数,值越接近于1越好；

Rmse表示均方根误差,值越小，表示拟合的曲线和真实数据越接近； 由结果可知sse和rmse两者的值都十分接近于0，即表示拟合的曲线和真实数据十分接近，并且确定系数和调整后的确定系数两者的值都接近于1，即综合表示拟合结果较好。

八、模型的评价与推广

本文运用质量守恒定律，利用微分方程模型，得到了低浓度一侧溶液中溶质的浓度与时间的微分方程模型，通过参数求解和简化得到了浓度与时间的指数关系。利用1stOpt软件和Matlab软件进行编程拟合，求得渗透率K。拟合结果可信度较高。但模型建立忽略了实验过程中A、B溶液体积的变化，这对最终结果的精确度有一定的影响。总的来说本模型简单、易于求解。

九、参考文献

[1] 姜启源.数学建模[M].北京：高等教育出版社，1991第二版. [2] 吴孟达，数学建模的理论与实践[M].国防科技大学出版社，1999.

十、附录

程序：

1stOpt软件程序：

Title \"分子穿透能力测试\"; Parameters a,b,k; Variable x,y;

Function y=a + b*exp(-0.2*k*x);

Data; 100 0.00454 200 0.00499 300 0.00535 400 0.00565 500 0.00590 600 0.00610 700 0.00626 800 0.00639 900 0.00650 1000 0.00659

Matlab程序：

clc; clear all;

t = [100:100:1000]; t = t';

C = 1e-03*[4.54,4.99,5.35,5.65,5.90,6.10,6.26,6.39,6.50,6.59]; C = C';

st = [0,0.0,0.01]; %初值

ft_ = fittype('a+b*exp(-0.2*K*t)', 'dependent', {'c'}, ... 'independent', {'t'}, 'coefficients', {'a', 'b', 'K'});

%f = fit(t,c,ft_ ,'Startpoint',st)

[f, goodness] = fit(t, C, ft_, 'Startpoint', st) figure

plot(f, 'predobs', 0.95); hold on

plot(t, C, 'b*')