【数学建模】九：MATLAB机理建模方法

在数学建模中，如果遇到一个非典型的数学建模问题（非数据、优化、连续、评价），那么这种情况下，通常需要用到机理建模方法了。

机理建模就是根据对现实对象特性的认识，分析其因果关系，找出反映内部机理的规则，然后建立规则的数学模型。机理建模的经典案例有很多，比如万有引力公式的推导过程。机理建模常见的有两类，一类是推导法机理建模，类似于微分方程建模，常用于动力学的建模过程，比如化学中反应动力学，还有各种场的方程，比如压力场、热场方程等；一类是包含一个或几个类别对象的复杂系统问题，常通过元胞自动机-仿真法来进行机理建模。下面将介绍这两类机理建模的具体 MATLB 实现过程。

1. 推导法机理建模

1.1 问题

某种医用薄膜有允许一种物质的分子穿透它（从高浓度的溶液向低浓度的溶液扩散）的功能，在试制时需测定薄膜被这种分子穿透的能力。测定方法如下：用面积为 S 的薄膜将容器分成体积分别为 VA、VB 的两部分，在两部分中分别注满该物质的两种不同浓度的溶液。此时该物质分子就会从高浓度溶液穿过薄膜向低浓度溶液中扩散。已知通过单位面积薄膜分子扩散的速度与膜两侧溶液的浓度差成正比，比例系数 K 表征了薄膜被该物质分子穿透的能力，称为渗透率。定时测量容器中薄膜某一侧的溶液浓度值，可以确定 K 的数值，试用数学建模的方法解决 K 值的求解问题。

1.2 假设和符号说明

为了便于建模，作以下几点假设：

① 薄膜两侧的溶液始终是均匀的，即在任何时刻膜两侧的每一处溶液的浓度都是相同的。

② 当两侧浓度不一致时，物质的分子穿透薄膜总是从高浓度溶液向低浓度溶液扩散。

③ 通过单位面积膜分子扩散的速度与膜两侧溶液的浓度差成正比。

④ 薄膜是双向同性的即物质从膜的任何一侧向另一侧渗透的性能是相同的。

图1 圆柱体容器被薄膜截面S阻隔

同时，约定需要用到的几个数学符号：

① CA(t)、CB(t) 表示 t 时刻膜两侧溶液的浓度;

② aA、aB 表示初始时刻两侧溶液的浓度 (单位:毫克/立方厘米);

③ K 表示渗透率;

④ VA、VB 表示由薄膜阻隔的容器两侧的体积;

1.3 模型的建立

考察时段 [ t , t+Δt ] 薄膜两侧容器中该物质质量的变化。以容器A侧为例，在该时段物质质量增加量：VACA(t+Δt)-VACA(t)

另一方面由渗透率的定义我们知道，从B侧渗透至A侧的该物质的质量为：SK(CB-CA)Δt

由质量守恒定律，两者应该相等，于是有：

VACA(t+Δt) - VACA(t) = SK(CB-CA)Δt

两边除以Δt，令Δt→0并整理得

dCA/dt = SK(CB-CA)/VA (3-1)

且注意到整个容器的溶液中含有该物质的质量应该不变，即有下式成立：

VACA(t) + VBCB(t) = VAaA + VBaB

CA(t) = aA + VBaB/VA - VBCB(t)/VA

代入式（3-1）得：

再利用初始条件 CB(0)=aB 解出：

至此，问题归结为利用CB在时刻tj的测量数据 Cj (j=1,2,…,N) 来辨识参数K和aA，aB对应的数学模型变为求函数：

问题转化为求函数

1.4 模型中参数的求解

例如，设VA = VB = 1000 立方厘米，S = 10 平方厘米，对容器的 B 部分溶液浓度的测试结果如表 1 所示。

表1 容器B部分溶液测试浓度

其中 Cj 的单位：毫克/立方厘米。

此时极小化的函数为：

下面用MATLAB进行参数求解：

1）编写M-文件(curvefun.m)

function f = curvefun(x,tdata)

f = x(1)+x(2)*exp(-0.02*x(3)*tdata)；%其中 x(1) = a;x(2) = b;x(3) = k;

2）编写程序(test1.m)

tdata = linspace(100,1000,10);

cdata = 1e-05.*[454 499 535 565 590 610 626 639 650 659];

x0 = [0.2,0.05,0.05];

opts = optimset( ‘lsqcurvefit’ );

opts = optimset( opts, ’PrecondBandWidth’, 0 )

x=lsqcurvefit (’curvefun’,x0,tdata,cdata,[],[],opts)

f= curvefun(x,tdata)

plot(tdata,cdata,’o’,tdata,f,‘r-‘)

xlabel(’时间(秒)’)

ylabel(’浓度（毫克/立方厘米）’)

3) 输出结果:

x =

0.0063 -0.0034 0.2542

% 即表示k = 0.2542, a = 0.0063, b =-0.0034时

f =

0.0043 0.0051 0.0056 0.0059 0.0061 0.0062 0.0062 0.0063 0.0063 0.0063

图2 模型拟合曲线与溶液实际测试浓度

曲线的拟合结果如图 2，进一步可求得：aB = 0.004（毫克/立方厘米），aA = 0.01（毫克/立方厘米）。

2. 元胞自动机-仿真法机理建模

元胞自动机 (Cellular Automata，简称 CA)，亦被称为细胞自动机。CA 的经典案例是定义一个网格，网格上的每个点代表一个有限数量的状态中的细胞。过渡规则同时应用到每一个细胞。典型的转换规则依赖于细胞和它的（4 个或 8 个）近邻的状态，虽然临近的细胞也同样使用。 CA 的应用在并行计算研究，物理模拟和生物模拟等领域。在数学建模中，一般是借鉴元胞自动机的概念，应用于具体的适合于机理建模的问题中。这类问题的典型特征是，所研究的问题是一个系统问题，系统是由若干个一个或几个不同类的对象组成，经典的模型不适应。典型的问题如滴滴打车问题（2015）, 开发小区问题（2016）。

这类问题，首先要分析系统内的对象，从微观角度研究每个对象的行为规则（模型），然后通过动态仿真研究系统内的对象随着时间或其他物理量的变化趋势，然后再根据目标综合评估系统。总结下来，实现步骤是：

1) 定义元胞的初始状态

2) 定义系统内元胞的变化规则

3) 设置仿真时间，输出仿真结果

对于这类仿真， MATLAB 的优势非常明显，下面通过介绍一个典型的 CA 的 MATLAB 实现过程：

%% 元胞自动机（CA）MATLAB实现程序

clc, clf,clear

%% 界面设计（环境的定义）

plotbutton=uicontrol(’style’,‘pushbutton’,…

’string’,‘Run’, …

‘fontsize’,12, …

‘position’,[100,400,50,20], …

’callback’, ‘run=1;’);

% 定义 stop button

erasebutton=uicontrol(’style’,‘pushbutton’,…

‘string’,‘Stop’, …

’fontsize’,12, …

’position’,[200,400,50,20], …

’callback’,‘freeze=1;’);

% 定义 Quit button

quitbutton=uicontrol(’style’,‘pushbutton’,…

‘string’,‘Quit’, …

‘fontsize’,12, …

‘position’,[300,400,50,20], …

‘callback’,’stop=1;close;’);

number = uicontrol(’style’,‘text’, …

’string’,‘1’, …

‘fontsize’,12, …

’position’,[20,400,50,20]);

%% 元胞自动机的设置

n=128;

z = zeros(n,n);

cells = z;

sum = z;

cells(n/2,.25*n:.75*n) = 1;

cells(.25*n:.75*n,n/2) = 1;

cells = (rand(n,n))<.5 ;

imh = image(cat(3,cells,z,z));

axis equal

axis tight

% 元胞索引更新的定义

x = 2:n-1;

y = 2:n-1;

% 元胞更新的规则定义

stop= 0; % wait for a quit button push

run = 0; % wait for a draw

freeze = 0; % wait for a freeze

while (stop==0)

if (run==1)

%nearest neighbor sum

sum(x,y) = cells(x,y-1) + cells(x,y+1) + …

cells(x-1, y) + cells(x+1,y) + …

cells(x-1,y-1) + cells(x-1,y+1) + …

cells(3:n,y-1) + cells(x+1,y+1);

% The CA rule

cells = (sum==3) | (sum==2 & cells);

% draw the new image

set(imh, ’cdata’, cat(3,cells,z,z) )

% update the step number diaplay

stepnumber = 1 + str2num(get(number,’string’));

set(number,’string’,num2str(stepnumber))

end

if (freeze==1)

run = 0;

freeze = 0;

end

drawnow %need this in the loop for controls to work

end

运行这段代码，可以得到如图 3 所示的初始图。

图3 元胞自动机初始图像

点击 Run 按钮，可以得到如图 4 所示的仿真图。

图4 元胞自动机仿真图像

如果改变运行规则，还可以到其他图像，如图 5 所示。

图5 元胞自动机仿真图像（更改规则后）

以上只是给出一个 MATLAB 实现典型元胞自动机的一个框架，具体建模的时候，还要根据具体问题，灵活定义元胞，更新规则，以及系统输出。比如在 CUMCM 2015 年打车问题中，元胞就是打车人和出租车，更新规则则是当打车人发出打车信号时，周边出租车的影响规则，系统输出则是评价指标。所以说元胞自动机只是一个概念，在实际建模问题中，还是要根据特定的问题再灵活运用。

【数学建模】九：MATLAB机理建模方法相关推荐

Matlab数学建模(九)：机理建模方法
一.学习目标 (1)掌握推导法机理建模 (2)掌握元胞自动机二.实战演练 1. 推导法机理建模 1.1 问题某种医用薄膜有允许一种物质的分子穿透它(从高浓度的溶液向低浓度的溶液扩散)的功能,在试制 ...
建模软件matlab,方程建模与MATLAB软件
第1章MATLAB的基本使用方法 1.1MATLAB概述 1.1.1MATLAB发展史 1.1.2MATLAB帮助 1.2MATLAB基础知识 1.2.1MATLAB变量 1.2.2MATLAB数据类 ...
数学建模常用方法 | matlab代码 | 二十三种数学建模方法 |2022赛前突击 |模型代码 |比赛比用、简单高效| 分享
为是赛前突击,所以就不过多的介绍理论知识了,直接上案例,matlab代码更加详细例题解析: 公众h:露露IT 目录 1.类比法 2.二分法 3.量纲分析法 4.图论法 5.差分法 6.变分法 7.数 ...
matlab建模总结,MATLAB 数学建模方法与实践（第 3 版）
本书从数学建模的角度介绍了 MATLAB 的应用,涵盖了绝大部分数学建模问题的 MATLAB 求解方法.全书共 5 篇.第一篇是基础篇,介绍基本概念,包括 MATLAB 在数学建模中的地位.数学模型的 ...
蹦极模型matlab仿真,科学网—蹦极的数学建模及其龙格-库塔法求解方法 - 赵也非的博文...
论文: 蹦极的数学建模及其龙格-库塔法求解方法在"华东师范大学首届研究生数学建模竞赛"中,获得二等奖. 发表日期: 2007年5月摘要: 本文通过参照题中给出的数据,对蹦极者在 ...
卓金武——从数学建模到MATLAB
卓金武--从数学建模到MATLAB 2013-9-4 09:48| 发布者: ilovematlab| 查看: 9647| 评论: 40 摘要: 人物简介--卓金武(Steven),MathWorks ...
最优化建模算法理论之Goldstein准则（数学原理及MATLAB实现）
文章目录一.前言二.Goldstein准则 1. 定义 2. 几何含义三.代码实现四.与Armjio准则的对比五.总结一.前言为了克服 Armijo 准则的缺陷,我们需要引入其他准则来保 ...
matlab app设计步骤_1.1数学建模与MATLAB–MATLAB入门
1.1数学建模与MATLAB–MATLAB入门关注本专栏,继续分享数学建模与MATLAB知识一.MATLAB是什么? MATLAB 是目前在国际上被广泛接受和使用的科学与工程计算软件.虽然 Cle ...
matlab或_数学建模与MATLAB——MATLAB入门
点击上方"蓝字",有更多精彩等着你噢! 关注本专栏,我们将继续分享数学建模与MATLAB知识. 你想要的,我都有! 一MATLAB是什么?MATLAB 是目前在国际上被广泛接受和使 ...

【数学建模】九：MATLAB机理建模方法

【数学建模】九：MATLAB机理建模方法相关推荐

最新文章

热门文章