（计算数学专业论文）药品扩散中的优化控制及其数值方法.pdf

d r u gd i f f u s i o ni no p t i m i z i n gc o n t r o la n di t sn u m e r i c a lm e t h o d b y x i a n gz i q u a n b e ( s h a o y a n gc o l l e g e ) 2 0 0 8 at h e s i ss u b m i t t e di np a r t i a ls a t i s f a c t i o no ft h e r e q u i r e m e n t sf o rt h ed e g r e eo f m a s t e ro fm a t h e m a t i c s c o m p u t a t i o n a lm a t h e m a t i c s c h a n g s h au n i v e r s i t yo fs c i e n c e t e c h n o l o g y s u p e r v i s o r p r o f e s s o rl iy o u y u n m a y ，2 0 1 1 长沙理工大学 学位论文原创性声明 本人郑重声明：所呈交的论文是本人在导师的指导下独立进行研究所取得的 研究成果除了文中特别加以标注引用的内容外，本论文不包含任何其他个人或 集体已经发表或撰写的成果作品对本文的研究做出重要贡献的个人和集体，均 已在文中以明确方式标明本人完全意识到本声明的法律后果由本人承担 作者签名：匈年松 日期矽，j 年夕月力夕日 学位论文版权使用授权书 本学位论文作者完全了解学校有关保留、使用学位论文的规定，同意学校保 留并向国家有关部门或机构送交论文的复印件和电子版，允许论文被查阅和借阅 本人授权长沙理工大学可以将本学位论文的全部或部分内容编入有关数据库进行 检索，可以采用影印、缩印或扫描等复制手段保存和汇编本学位论文同时授权 中国科学技术信息研究所将本论文收录到中国学位论文全文数据库，并通过网 络向社会公众提供信息服务 本学位论文属于 1 、保密口，在年解密后适用本授权书 2 、不保密曲 ( 请在以上相应方框内打“”) 作者签名：么勿 红日期：刀卅年岁月加 导师签名：多毒饭弓 日期： l f 年乡月z 日 摘要 许多实际问题中的物理、医学、化学、生物过程都可以用扩散方程来描述 为了证明药物有效成分的溶解与扩散对药物的生物利用度以及药效和研究药物载 体对药物溶解与扩散性质的影响，从而需要确定合适的扩散系数，本文就是以此 为口的展开研究的针对药品扩散问题中的扩散系数难于获取的问题以及溶液本 身及药品可溶性，本文提出了对不同容器形状的优化控制模型及其相应的数学方 法来反演模拟过程，从而提高了药物的安全性、实用性和有效性 本文具体结构如下： 第一章为绪论，主要介绍了本论文的研究背景和选题依据，以及研究内容、 方法和创新 第二章提出了一种估计圆柱体传递装置模型的有效扩散系数的优化方法和一 种估计圆柱体传递装置中药物传递模型的依赖时间的有效扩散系数参数的优化方 法实例都说明了算法的有效性 第三章针对二维圆盘传递装置和球形传递装置的线性扩散问题，利用分离变 量方法和最小二乘法，得到了线性扩散的解析解，实例说明了算法的收敛性和有 效性 第四章针对二维圆盘传递装置的非线性扩散和球形传递装置的非线性扩散问 题，利用线性扩散去拟合非线性扩散，得到了非线性扩散的半解析解，结果显示 了算法的有效性 第五章总结药物线性扩散和药物非线性扩散方法的优点和缺点，并展望下一 步工作 关键词：有效扩散系数；药物控制释放；g a u s s n e w t o n 方法；扩散系数；最优 化算法 a b s t r a c t m a n yp r a c t i c a lp r o b l e m si np h y s i c s ，m e d i c i n e ，c h e m i c a l ，b i o l o g i c a lp r o c e s sc a nb e u s e dt ob ed e s c r i b eb yt h ed i f f u s i o ne q u a t i o n i no r d e rt op r o v et h a tt h ed r u ge f f e c t i v e i n g r e d i e n t so fd r u g sa n ds p r e a d t h eb i o a v a i l a b i l i t ya n de f f i c a c yo fd r u gr e s e a r c h ， d e t e r m i n et h ea p p r o p r i a t ei n v e r s i o nd i f f u s i o nc o e f f i c i e n t ，t h i sp a p e ri sd o n e b e c a u s e i ti sd i f f i c u l tt oo b t a i nt h ed i f f u s i o nc o e f f i c i e n to ft h ed r u gd i f f u s i o np r o b l e m ，w es h a l l p r o p o s e t h e o p t i m i z a t i o n c o n t r o lm o d e li nt h ed i f f e r e n tc o n t a i n e r s a n dt h e c o r r e s p o n d i n gm a t h e m a t i c a lm e t h o d sf o ri n v e r s i n gt h es i m u l a t i o np r o c e s s t h es t r u c t u r eo ft h i sp a p e ri so u t l i n e da sf o l l o w s ： t h ef i r s tc h a p t e ri sa ni n t r o d u c t i o n ，w h i c hd e s c r i b e st h er e s e a r c hb a c k g r o u n d ，t h e m o t i v a t i o no ft h e s i sa sw e lla st h er e s e a r c hc o n t e n t ，m e t h o da n di n n o v a t i o n t h es e c o n dc h a p t e rp r o p o s e st w oo p t i m i z a t i o nm e t h o d s o n ei st oe s t i m a t et h e d i f f u s i o np a r a m e t e r si nt h ec y l i n d e rt r a n s f e rd e v i c em o d e l ，t h eo t h e ri st oe s t i m a t et h e e f f e c t i v et i m ed e p e n d e n td i f f u s i o nc o e f f i c i e n tp a r a m e t e ri n t h ec y l i n d e rt r a n s f e r d e v i c e e x a m p l e ss h o wt h a tt h ep r o p o s e dm e t h o d s a r ev a l i d a i m e da tt w od i m e n s i o n a ld e v i c e s ：t h ed i s ca n ds p h e r i c a l ，s o m em e t h o d sb a s e d0 1 1 t h es e p a r a t i o nv a r i a b l e si d e aa n dt h el e a s ts q u a r em e t h o di nt h et h i r dc h a p t e ra r eg i v e n t oc o m p u t et h el i n e a rd i f f u s i o np r o b l e m s w eo b t a i nt h ea n a l y t i c a ls o l u t i o n so ft h e l i n e a rd i f f u s i o np r o b l e ma n ds o m ee x a m p l e si l l u s t r a t et h ei rc o n v e r g e n c ea n dv a l i d i t y i nt h ef o u r t hc h a p t e rf o rt h en o n l i n e a rd i f f u s i o np r o b l e m si nt w od i m e n s i o n a l r e l e a s i n gd i s c o rs p h e r i c a ld e v i c e ，s u b s t i t u t i n gt h el i n e a rd i f f u s i o np r o c e s sf o rt h e n o n l i n e a rd i f f u s i o np r o c e s s ，w eo b t a i nt h es e m i - a n a l y t i c a ls o l u t i o no ft h en o n l i n e a r d i f f u s i o np r o b l e m t h er e s u l t sa l s os h o wt h ee f f e c t i v e n e s so ft h ep r o p o s e da l g o r i t h m i nt h ef i f t hc h a p t e rw es u m m a r i z et h ea d v a n t a g e sa n dd i s a d v a n t a g e so ft h e m e t h o d sf o rt h el i n e a rd i f f u s i o np r o b l e ma n dn o n l i n e a rd i f f u s i o np r o b l e m so ft h ed r u g r e l e a s i n gp r o c e s s a tl a s t ，s o m ef u t u r ew o r ki sl i s t e d k e yw o r d s ：e f f e c t i v ed i f f u s i o nc o e f f i c i e n t ：d r u gr e l e a s i n g ；g a u s s - n e w t o n m e t h o d ；d i f f u s i o nc o e f f i c i e n tso p t i m i z a t i o na l g o r i t h m 目录 摘要i a b s t r a c t i i 第一章绪论 1 1 本文研究的背景( 1 ) 1 2 国内外研究现状( 1 ) 1 3 本文研究的意义及方法( 2 ) 第二章三维圆柱体传递装置的药物线性扩散的有限差分方法 2 1 与时间无关的优化控制模型及算法( 4 ) 2 2 依赖时间的优化控制模型及算法( 1 5 ) 第三章药物线性扩散的解析方法 3 1二维圆盘传递装置的药物线性扩散的优化控制模型及算法( 3 0 ) 3 2 球形传递装置的药物线性扩散的优化控制模型及算法( 4 1 ) 第四章药物非线性扩散的半解析方法 4 1二维圆盘传递装簧的药物非线性扩散的优化控制模型及算法( 5 1 ) 4 2 球形传递装置的药物非线性扩散的优化控制模型及算法( 6 1 ) 第五章结果分析与结论( 6 9 ) 参考文献( 7 0 ) 致谢( 7 3 ) 附录( 攻读学位期间发表论文目录及参与科研项目) ( 7 4 ) 1 1 本文研究的背景 第一章绪论弟一早瑁t 匕 物理、医学、化学和生物领域中的许备模型都可归结为所谓的关于反演扩散 系数的反应扩散方程现代科学技术的发展很大程度上是南于物理、化学和生物学 的成就和发展，并且这些科学自身的精确化又是他们取得发展的重要保证科学的 精确化往往是通过建立数学模型来实现的近年来，药物扩散方程研究受到广泛重 视，科技研究论文层出不穷这是由于药物扩散方程具有强烈的应用背景在数学 理论研究巾，扩散方程解可以揭示方程本身许多重要的性质在实际应用巾，本 文研究在一个容器中，装有某种溶液，把药剂放进去，就可以模拟药剂在人体中 的吸收情况药剂扩散慢，则达不到治疗的效果；药剂扩散快，则对人体有害因 此，药物有效成分的溶解与扩散对药物的生物利用度以及药效有重要的影响基 于此，本文建立了一个研究药物扩散的反演扩散系数刀的线性模型( 1 ) 及非线性 模型( 2 ) 掣一刚 可0 2 c + 石o c 歹i + 万0 2 c 7 1 + 害 - 0 ， 吣只力q 皇9 2 堡型：o o r 。 叫厂，只石o ) = h ( ，) ， ， o ，( 乃绣刁a q ( 1 ) ( ， 绣才q 掣一v ( 乃( q v 印= 0 ， ， o ，( ，只力q 坌三( 三堡墨垒：o ， o r “以互o ) = h ( 厂) ， t o ，( 层刁a q ( 2 ) ( 层力q 其巾刀是扩散系数，f 是未知浓度，q 、q 。是求解区域 1 2国内外研究现状 针对药品扩散问题中的扩散系数难于获取的问题，本文提出了针对不同容器 形状的优化控制模型及其相应的数值方法，现对国内外研究现状给予叙述 程爱杰等以多孔介质中的控制释放由边界积分一常微分方程捕述，溶质迁移 由带第三类边界条件的对流扩散( 含机械弥散) 方程描述，构造了这一非线性耦合 问题的有限元半离散格式，利用先验估计理论进行了收敛性分析 兰晓林乜1 对一类非线性反应扩散对流方程的有限差分方程组给出了一类数值 计算方法，包括时问依赖问题的有限差分方程组和相应的定常问题的有限差分方 程组 芦碧波口3 主要研究了基于偏微分方程( p d e ) 的处理方法，特别是非线性扩散 模型在图像去噪中的模型和图像分割中的应用 黄国灿h 1 认为物理，化学和生物领域中的许多模型都可归结为所谓的反应扩 散方程反应扩散方程有一类重要的解，就是形如k ，) = p + 卅的解在数学理论 研究中，行波解可以揭示方程本身许多重要的性质；在实际应用中，行波解可以 很好的解释自然界中的一些波的传播 p a v e l 晦6 1 等人研究了小分子在无定形聚合物中气体的扩散系数，验证了自 南体积在扩散系数巾的重要性b o y d n3 等人则模拟了阻挡材料巾气体的扩散系数， 发现与其他同类材料相比，阻挡材料中气体的扩散系数较小h o f m a n n 睁1 们等人研 究了小分子在玻璃态和高弹态聚合物中扩散的不同m u l l e r l e r p l a t h e n 等人也 研究了气体在聚合物中的扩散v a nd e rv e g t n 2 1 讨论了温度对于二氧化碳在聚合 物中扩散的影响s u t e r 1 3 1 4 1 采用不同的方法，模拟聚合物中气体的扩散过程l u n 司 等人研究了亲水性链在疏水表而的行为g r e e n f i e l d u 町等人用r k m c 方法，研究聚 合物中的扩散行为y a r o u s k y n 7 3 等人研究了交联性聚合物中的分子扩散行为 w a n g 等人n8 1 钔研究发现，聚合物巾大的自南体积有利于小分子的扩散目前，很 多其他的课题组也进行了这方面的研究工作 2 0 - 4 4 1 3 本文研究的意义及方法 扩散和对流扩散过程出现在许多领域，比方说地球地理学、工程学、经济学、 金融学，生物医学科学等 4 5 - 4 8 在许多情况下，扩散系数未知，需要我们根据实验 或者寻找探索性数据来确定这个扩散系数虽然扩散过程由函数空间、时问及物 质浓度控制，但是我们希望可以找到一个适合的扩散系数来满足这个过程 2 在过去的几十年，药物控制装置设计已经吸引了很多注意力，装置的功能和 性能对装置的有效扩散率是很关键的虽然药物装置系统的扩散系数主要由孔隙 率和材料的其他一些性质决定，当这些性质已知，如何选取有效扩散系数成为一 个主要问题有各种各样的现有技术来确定有效扩散系数的阳离子x l o u ， s w a n g 。等以为这些技术依据经验或半经验模型，这些模型依据药物传递机制或 者二维或特殊情况的扩散方程的解析解 4 9 0 s o 在实际应用中，三维模型装置经常用 到如果仅仅依靠经验或半经验模型，选择扩散系数是一件很困难的事情然而， 对于解析解，仅限于二维模型或特殊装置晡1 哪! 在文献嘲1 的基础上，为了更好地分析三维模型，本文提出一个数值方法去获 取扩散和对流扩散过程的扩散系数然而，因为不断扩散，很难获得最佳逼近的质 量比，冈为药物的释放速率在不同时间不同过程是不同的由于这些原冈，在这篇 论文中，针对扩散过程和非线性扩散过程，我们将提出一种新的数值方法去获取 扩散系数通过研究药物在不同材料中的扩散性质，可以提高药物的安全性、有效 性和适应性因此，可以指导人们对不同药物选取合适的控释材料，因此有非常重 要的现实意义 在本课题研究中，重点是阅读大量文献，以及了解在这个领域巾前人所作的 成果及研究方法本文的结构如下：一个是药品扩散问题中的反演扩散系数的优化 控制线性模型及其算法另一个是药品扩散问题中的反演扩散系数的优化控制非 线性模型及其算法下一章提出了一种估计圆柱体传递装置模型的有效扩散系数 的优化方法及一种估计圆柱体传递装置中药物传递模型的依赖时间的有效扩散系 数参数的优化方法实例都说明了算法的的有效性第三章针对二维圆盘传递装 置的线性扩散和球形传递装置的线性扩散问题，利用分离变量方法和最小二乘法， 得到了线性扩散的解析解，实例说明了算法的收敛性第四章针对二维圆盘传递装 置的非线性扩散和球形传递装置的非线性扩散问题，利用线性扩散去拟合非线性 扩散，得到了非线性扩散的半解析解，结果显示了算法的有效性第五章总结了药 物线性扩散和药物非线性扩散的方法的优点和缺点，并展望下一步工作 3 第二章三维圆柱体传递装置的药物线性扩散的有限 差分方法 本章主要讨论与时间无关的反演扩散系数和依赖时问的药物扩散系数的优化 控制线性模型及算法此章共分两节 2 1 与时间无关的优化控制线性模型及算法 2 1 1 模型介绍 针对三维情形，装置为圆柱体容器，装置由大、小两个圆柱体容器构成，小 圆柱形容器同定在大圆柱形容器内部，其中小容器外面和大容器里面盛满了液体， 装置见图1 小圆柱形容器和大圆柱形容器的半径、高、包含的区域分别表示为百， 属，f l 。和呢，鹰，f l ，其中小圆柱形容器内部有数量为m o 的药物考虑如下扩散 系数为参数的筲卡尔坐标扩散方程： 皇掣一d s c ( 五另弓力：0 8 t “；1 。 型些y - 0 8 n ； “五另互o ) = h ( 石另力， 其中，d 表示常数：q 五月石，) 表示未知药物浓度 ， o ，( 石另习q ， o ，( 五月刁a q ( 工男刁q l 图l ：圆柱体容器传递装置 4 ( 2 1 1 ) 关于初始条件h ( 五另刁，本文假定在，= 0 时刻，药物浓度一致并且液体巾的浓度 为0 ，即 毗胁黪舢如腓q ( 2 1 2 ) 1 0( 五男力q 、q l 因此，扩散过程确定了这个扩散系数，换言之，等价解决下列的优化问题 问题1 寻求刀满足 嬲 ( m 。- m ? ) 2 + ( m ：一m ：) 2 + - + ( m ，一m ：) 2 ) ， ( 2 1 3 ) 其中，m ? ，m ：，m ：表示给出的实验数据，m 。，m ：，m ，由下列的计算公式给 出： m ，= ，h m 。q 五男磊) 翻掀，= l ，2 ，弓 ( 2 l 4 ) 这里，“五另互，) 表示方程( 2 1 1 ) 的解 为了解出c ，本文在极坐标中求解问题1 因此，问题l 就等价转化为下列 的问题2 问题2 寻求刃满足 lim(m胁o。一m ? ) 2 + ( m ：一m ：) 2 + 十( m ，一m ：) 2 ( 2 1 5 ) 其中，m 。，m ：，m ，由下列的计算公式给出： m ，= 觚、na ( r , 0 , z , 1 ) d r d o d z , i = 1 ，玑岛 ( 2 1 6 ) 并且“，只互，) 由下列的极坐标扩散方程确定： 垫o 型t 一刃黟+ 等圭+ 等专+ 翻= o ， 0 ，( 删q 皇掣：0 ， o ，( ， 层力a q ( 2 1 7 ) a ， 。 “ 只五o ) = h ( 只刁( 只力q h c ，只力= 手之乏三三三三毫i 鹰 5 ( 2 1 8 ) 2 1 2离散方程 为了获取最佳扩散系数刃，本节采用一个有限差分方法去离散并求解偏微分 方程( 2 1 7 ) 在极坐标扩散方程( 2 1 7 ) 中，设汐，z 的范围分别表示为( 0 ，2 力，( 0 ，呓) ， ( 0 ，甸，它们划分的相等子区问的个数分别为p ，口r ，相应步长分别为a o , a t , 止， 其中p ，口r 都为正整数因此，就在装置内部定义了一个空问网格区域，网格区 域由网格点组成，即由纪，乃，乞组成，其中，= l ，p ，一，= l ，口七= l ，刀给定 一个时间步长及一个正整数t ，当= 1 ，2 ，t 时，令= 卜1 ) a t 利用时间和 空间的离散数据，本节在偏微分方程( 2 1 7 ) 中得到了下列偏导数的近似计算： 一o c 二二 ( 2 1 9 ) 一一 l 二上了， 丝互业二曼占 8 ra r 等丛呜笋塾a ， 等盈嚆笋坠8 伊铲 等亟嚆笋鱼8 方时 ( 2 1 1 0 ) 将式( 2 1 9 ) 一( 2 1 1 ：3 ) 代入方程( 2 1 7 ) ，整理可得偏微分方程组： 一日4 t 一4 吃4 一4 _ + 忍勃一忍吼_ t 一忍乞甜一匀乞t m2 警 d , , q - l ，t 一岔t 2 0 吒鼻一c ：劓= o ( 2 1 1 4 ) 上。一, i j - - 0 c 幺t = h ( ，只刁 对于所有的z 磊都成立其中4 ，皿，4 表示为( 2 1 1 5 ) 很显然，关于乞。的线性方程组对于所有的z 毛都满足 在上述离散方法中，本文虽然采用均匀分割，简洁符号来表示，很显然，用 不均匀步长分割的离散方法也是适合的 6 4 = 面d 露， 4 = 参， 。 忍= 参， 忍= 。侈+ 南+ 南击b 1 ， 忍= 蕊d ， 色= 侈+ 毒】， 4 = 告 2 1 3 解决最小二乘法问题的数值方法 本节将提出关于问题2 的一些算法令 e ( 功= ( m 。一m r ) 2 + ( m ：一a f ：) 2 + + ( i ，一h i ：) 2 三( 肜一m ) t 彳( 肜一m ) ( 2 1 1 6 ) 其中，m = ( h i 。( 切，m ：( 切，m ，( 功) t 及盯= ( a i ：，m 弘，a i ：) t 考虑问题2 的数值方法，从初始估计开始，采用迭代方法解决问题2 其 每一步的增量为，使得e ( + ) 最小，其中和分别表示刀的第，步近 似值和第，步增加值当l l m m l i i m 0 小到一个给定的正常数时，迭代过程终 ，一 止 为了计算每一步的增量，根据文献 1 4 巾的思想，建立了一个 g a u s s n e w t o n 方法 向量值函数的泰勒公式为 m = m 7 + 矽+ 兰 ( ) 2 刁， 其中 z 警，訾，訾】t 7 g 表示m 在点十肛处的二次导数向量值，其中0 _ s a si ，本节取夕= i 根 式( 2 1 1 7 ) ，去掉二次项得到 m = m i 七j 6 d 其中，当很小时，e ( + 桫) 可以近似表示为 e ( + 尻) = ( 1 7 + z 五莎一m ) t a ( m 7 + z 五一m ) = ( a 1 7 一m 。) t a ( m7 一m ) + ( l7 一m + ) t a z b l ( 2 1 1 + ( 彬) a ( m 7 一m ) + ( 矽) t a ( 删) 这是一个关于# l y 的二次函数，令二次函数取极值的最小点为尻矿，则二次函数满 足 v e ( + ) = o 其可以导出 ( z ) t a ( z ) 万= 一z w a ( 1 w - m ) ( 2 1 2 0 ) 式( 2 1 2 0 ) 的解表示第，步的搜索方向，即g a u s s - n e w t o n 方向根据方程 ( 2 1 2 0 ) 可以导出 = 一( ( z ) ta ( 硝1 a ( m 7 一m ) ( 2 1 2 1 ) 然后下一步的扩散系数d 可近似表示为 l + 1 = + 历矿( 2 1 2 2 ) 根据方程( 2 1 2 1 ) ，只要算出所有的z 关于的偏导数，就能获取第，步的 搜索方向根据连贯性，本节提出了一个计算，的算法 根据方程( 2 1 6 ) ，得到下列方程： 等= 几焉掣一，巴 汜m 3 ， 根据方程( 2 1 2 3 ) ，只有解出皇墨幽，才能求出皇坚：f 堡垒至墨望对于计 a d o d 。 算的每一个扩散系数，令乃= ( 1 + 功，其中万是一个很小的常数根据扩 散方程( 2 1 1 ) 和第，步步长的扩散系数及d ，就能获取的药物浓度c 及 勿的药物浓度e 利用下列偏导数公式近似计算： 皇g 堡2 墨墨型呈蔓， o ddx 6 7 8 其可以近似代替偏导数兰掣 可由下列公式算出： 因此，根据方程( 2 1 2 3 ) ，o c ( _ o , r 菘, z - , t , z r ) 等= f 7 上丘e 锸一川囊，只 可以看出，利用前面的有限差分方法可以解出偏导数垒学如果可 以求解偏导数璺墨掣的值，就可以根据方程( 2 1 2 3 ) 近似取偏导数 8 d 皇g 堡三墨垒望的值 根据前而的数值方法，下而给出求解最佳扩散系数刀的算法 在，= l ，2 ，3 ，f 时，设m ，表示相应时间的数值方程的质量点，m ：表示相 应时间的实验测量质量点根据下列的最4 , - 乘法算法( 厶吻，可以获取最佳扩 散系数刃 算法f l s a ) 1 选择一个正整数，令e 删为一个常数，让，= l ，j = 1 ，给出刀的初始值，转 第2 步 2 如果p ，根据式( 2 1 1 4 ) 求解扩散方程，以获取时间间隔【o 。，】的药物 浓度，转第3 步，否则，转第5 步 3 计算时间间隔i o ，l 的药物浓度导数，转第4 步 4 根据式( 2 1 6 ) ，利用扩散系数计算m ，令j = j + l ，转第2 步 5 根据方程( 2 1 2 3 ) ，计算扩散质量关于扩散系数的导数，转第6 步 6 根据式( 2 1 1 6 ) 计算e ( 功，如果e ( 功 ，转第9 步，否则，转第7 步 7 计算偏导向量霎当，形成偏导向量z ，根据式( 2 1 2 1 ) ，得到第，步的搜索 8 d 1 方向，转第8 步 8 令力= 刀+ 历7 ，= 件l 及= l ，转第2 步 9 输出摄佳扩散系数刀，算法终l 七 2 1 4 数值实例 为了验证数值方法的有效性，下面给出一个数值试验，在数值模拟试验中， 采用了不均匀步长分割的离散方法 9 测试问题如图1 所示，小圆柱形容器固定在大圆柱形容器内部，它们的半径和高见 表格i 表格2 表示的是药物释放传递装置中不同时问点的实验研究释放数据m ：m 本文为了节约计算时间，采用初始点d o = 1 1 0 _ 7 ，网格参数，= 1 0 0 x ，厂= 0 0 0 2 c m 及 a o = 0 8 0 n 解决此优化问题表格3 到表格8 表示的是前6 次迭代的结果然后，减少 网格尺寸，网格参数为a t = 5 0 s ，= 0 0 0 1 c m ，业= o 0 2 c m ，a o = 0 8 c m ，利用第6 次迭代结果么= 2 5 2 7 6 x 1 0 及初始估计继续迭代另外4 次表格9 到表格1 l 表示的是 m ，m 。的迭代结果 数值实验中，令每一次的总误差为 1 2 总误差( ，) 2 ( 计算数据一实验数据( 力) 2 户l 表格1 2 、表格1 3 分别表示每次的距离和每次的总误差图2 和图3 分别表示计 算的扩散系数和释放曲线根据表格1 2 、表格1 3 、图2 及图3 ，可以看出基于f d 的最优算法收敛从表格1 3 的最后一行可以得到，最小二乘法的误差最小图4 绘制了拟合曲线，表示计算结果m ，m 。非常接近实验数据m ：m 。因此，通过 对计算数据的分析，基于f d 的最优算法是有效的 表格1 ：药物释放传递装置尺寸 容器种类半径高 容器种类 2 8 3 9 9 c m1 0 1 9 c m 小圆柱形0 4 8 0 0 c m1 0 1 9 c m 表格2 ： 药物释放传递装置中的不同时间点的实验释放数据m ：m 。 时间( 秒) f h i 时间( 秒 m ! m 时间( 秒 m ：a 互 l8 0 00 1 9 7 4 7 01 0 8 0 00 3 9 3 5 3 61 1 7 7 2 00 7 7 9 3 0 5 3 6 0 0 0 2 4 2 9 8 8 1 6 2 0 0 0 4 4 4 0 6 6 1 8 3 4 2 0 0 8 8 2 8 6 3 5 4 0 00 3 4 2 1 5 72 4 7 2 00 5 0 6 9 8 8 2 0 0 8 2 0 0 9 0 9 4 0 5 7 2 0 00 37 0 6 4 28 9 5 2 00 7 5 1 2 6 5 2 6 2 0 8 0 1 0 0 0 0 0 0 1 0 表格3 ：药物释放传递装置中的不同时间点的第1 次迭代优化算法的释放数据 m ，m 。 ， l2 3 4 56 o 0 1 7 50 0 2 8 2o 0 3 7 20 0 4 5 2o 0 5 8 90 0 7 6 l 7891 01 11 2 0 0 9 8 30 2 0 4 4o 2 3 7 00 2 9 9 lo 3 1 3 30 3 5 8 4 表格4 ：药物释放传递装置中的不同时间点的第2 次迭代优化算法的释放数据 m ，m 。 12 3 4 56 0 0 5 4 60 0 8 4 70 1 0 8 10 1 2 7 80 1 6 1 00 2 0 1 5 7891 0l l1 2 0 2 5 2 90 4 8 0 1 0 5 4 2 00 6 4 5 80 6 6 6 90 7 2 6 7 表格5 ：药物释放传递装置中的不同时间点的第3 次迭代优化算法的释放数据 m ，u 。 l2345 6 0 0 9 4 70 1 4 2 1o 1 7 8 50 2 0 9 00 2 5 9 70 3 2 0 6 7891 01 l1 2 0 3 9 5 70 6 8 3 00 7 4 3 40 8 2 9 70 8 4 5 30 8 8 6 5 表格6 ：药物释放传递装置中的不同时间点的第4 次迭代优化算法的释放数据鼍 123456 0 1 2 5 30 1 8 5 60 2 3 1 50 2 6 9 70 3 3 2 50 4 0 6 6 7891 01 11 2 0 4 9 5 30 7 8 5 20 8 3 5 70 9 0 3 i0 9 1 4 70 9 4 4 8 表格7 ：药物释放传递装置中的不同时问点的第5 次迭代优化算法的释放数据篑 12345 6 0 1 3 7 60 2 0 2 90 2 5 2 50 2 9 3 60 3 6 0 80 4 3 9 4 7 8 91 0l l 1 2 0 5 3 2 0o 8 1 5 80 8 6 2 30 9 2 3 2 0 9 3 3 60 9 6 0 3 表格8 ：药物释放传递装置中的不同时问的第6 次迭代优化算法的释放数据鼍 l23456 0 i 3 0 00 1 9 1 50 2 3 8 10 2 7 6 80 3 3 9 90 4 1 3 7 7891 0l l1 2 0 5 0 0 50 7 6 3 30 8 0 6 0o 8 6 1 9o 8 7 1 40 8 9 6 0 表格9 药物释放传递装置中的不同时间点的第7 次迭代优化算法的释放数据卺 l23456 0 1 4 8 80 2 1 7 90 2 6 9 9o 3 1 2 80 3 8 2 30 4 6 2 0 7891 01 11 2 0 5 5 2 50 8 0 1 00 8 3 8 10 8 8 5 70 8 9 3 7o 9 1 4 1 表格1 0 ：药物释放传递装置中的不同时问点的第8 次迭代优化算法的释放数据 m ， m 。 l2 3 456 0 1 5 1 50 2 2 1 70 2 7 4 40 3 1 7 90 3 8 8 20 4 6 8 7 7891 01 11 2 0 5 5 9 50 8 0 5 70 8 4 2 10 8 8 8 60 8 9 6 40 9 1 6 3 1 2 表格1 1 ：药物释放传递装簧中的不同时间点的第9 次迭代优化算法的释放数据 m ， m 。 l 2 3456 o 1 5 1 60 2 2 1 80 2 7 4 50 3 1 8 00 3 8 8 30 4 6 8 8 7891 01 l1 2 0 5 5 9 70 8 0 5 80 8 4 2 20 8 8 8 70 8 9 6 50 9 1 6 3 表格1 2 ：迭代步长距离及药物释放传递装置巾的优化算法误差 迭代次数 l23 4 56 6 d4 3 8 8x 1 0 。7 7 5 9 l 1 0 。77 9 7 2 x1 0 73 7 0 7 x 1 0 4 6 4 1x 1 0 8 2 4 2 9 3x1 0 9 d0 5 3 9 l o 1 2 9 8 1 0 - 62 0 9 5 x1 0 42 4 6 6 x1 0 - 62 5 3 0 1 0 - 62 5 2 7 6x1 0 。6 误差 2 3 9 2 80 6 5 6 2o 1 4 5 00 0 4 4 00 0 3 5 7o 0 3 5 6 表格1 3 ：迭代步长距离及药物释放传递装置r f l 的小尺寸优化算法误差 迭代次数 7891 0 6 d 6 8 5 0 4 1 0 7 1 0 4 8 1 0 。71 7 6 7 8 x 1 0 91 3 2 6 6 x1 0 1 0 d 3 2 1 2 0 1 0 巧3 3 1 7 1 0 6 3 3 1 8 8 1 0 。6 3 3 1 9 0 x 1 0 。6 误差 0 0 4 4 20 0 2 7 70 0 2 7 40 0 2 7 3 图2 ：基于肋优化算法，从圆柱体药物容器巾获取的药品扩散系数 图3 ：基于f d 优化算法，每个迭代过程的步长距离 图4 ：基于肋优化算法，第1 。次迭代的实验数据翥苎和计算数据害三的比较 1 4 2 2 依赖时间的优化控制模型及算法 2 2 1 模型介绍 针对三维情形，装置为圆柱体容器，装箕由大、小两个圆柱体容器构成，小 圆柱体容器同定在大圆柱体容器内部，其中小容器外面和大容器里面盛满了液体， 装置见图1 小圆柱体容器和大圆柱体容器的半径、高、包含的区域分别表示为百， 属，q 。和五，鹰，q ，其中小圆柱形容器内部有数量为m o 的药物考虑如下的依 赖时问的扩散系数参数的极坐标扩散方程 掣叫“等+ o _ c l + 万0 2 c 7 1 + 期- o ， o ，( 删q o c ( , - , a , - - , o ：0 ， o ，( ，绣刁o f ( 2 2 1 ) or 。、。7 “只五o ) = h ( 厂) ， ( 只习q ， 其中，刀( ，) 表示时问的一个函数；q 五月石力表示未知药物浓度 关于初始条件，本文假定在，= 0 时刻，药物浓度一致并且液体中的浓度为0 ， 即 h ( ，力- j 等。o ( 棚刎q ( 2 2 2 ) 【0( r , o ，石，) 峨 凶此，寻求扩散参数刃( ，) 就等价于解决下列的最优化问题 问题3 寻求d ( ，) 满足 m 棚 ( m 。( ，) 一m ? ) 2 + ( m ：( ，) 一m ：) 2 + + m ，( ，) 一m ：) 2 ( 2 2 3 ) 其中，m 。( 力，m ：( 力，m ，( ，) 由下列的计算公式给出： m m = fr 。r q 以乞力搠幺川，2 ，3 ，p ( 2 2 4 ) 并且c ( r , o ，石，) 由极坐标扩散方程( 2 2 i ) 确定 2 2 2 采用有限差分公式求解极坐标扩散方程 为了获取依赖时间的最佳扩散系数参数d ，本节采用一个有限差分方法去离 散并求解偏微分方程( 2 2 1 ) 在极坐标扩散方程( 2 2 1 ) 巾，设厶，h = l ，2 ，t 表示扩散系数的离散次数。 i = i ，2 ，t 表示时间的离散次数，i = i ，2 ，p 表示0 的离散次数，= 1 , 2 ，夕表示 半径的离散次数，后= 1 ，2 ，r 表示坐标z 方向的离散次数，觚8 ，止分别表示 时间，、半径，、弧度日和纵坐标z 方向的步长利用时间和空间的这些离散数据， 本节在偏微分方程( 2 2 1 ) 中得到了下列偏导数的近似计算： a c d - c - 1( 2 2 5 ) 一一 二厶d ， 8 t 丝鱼型二鱼 a r 堡堡丛：! 笔：兰二三竺兰 8 ， 丝垒! ：丛竺! ：生生二三星丛 a 口2a 0 2 堡g 竺丝竺望二兰堡丛 a 心 将式( 2 2 5 ) 一( 2 2 9 ) 代入方程( 2 2 1 ) ，整理可得偏微分方程组： 一日1 k - - 忍h 一日 “+ 忍一忍矗。t 一忍 一日c 。纠2 号笋 ( 口一c b = 0 c i 伸一c 。? m = 0 c 印一乞l = 0 艺广h ( r , e ，刁 其中4 ，4 ，4 表示为 1 6 ( 2 2 6 ) ( 2 2 7 ) ( 2 2 8 ) ( 2 2 9 ) 忍= 黠， 4 = 辔， 4 = 辔， 口叫厶) 【吉+ 专+ 南+ 剖嘧l 忍= 黠， 色 辔+ 剖， 4 = 辔 2 2 3 从扩散过程中获取扩散系数参数的最佳方法 本币将给出解决同题3 的优化万法令 e ( 功= m 。( ，) 一m o ) 2 + ( m ：一m o ) 2 + + ( m ，( 力一m ：) 2 ， ( 2 2 1 2 ) 如果令彳为一个一阶单位矩阵，e ( 功就可转化为下列方程： e ( 功= ( 叫，) 一m ) t 彳( 川力一m ) ( 2 2 1 3 ) 其中，m = ( m ，( d ( ，) ) ，m ：( 刃( ，) ) ，m ，( 刀( ) ) ) t 及m = ( m ：，m ：，m ：) t 考虑问题3 的数值方法从初始估计= ( ( ) ，心) ，( o 。) ) t 开始，运 用迭代方法解决问题3 其每一步的增量为占= ( 6 g ) ，6 g ) ，6 ( 钿) ) ， 使得e ( + 6 ) 最小，其中和彩分别表示刃的第，步近