MATLAB在化学中的应用

硕士研究生课程论文

课程名称：MATLAB及其应用

题目：MATLAB在热物理学中的应用学院：

专业名称：

姓名：

学号：

任课教师：

提交时间：2013 年 11 月 21 日

MATLAB在热物理学中的应用

摘要：本文阐述了基于MATLAB的数值计算、可视化图形处理、开放式以及可扩充体系结构的特点,并介绍了高性能语言 MATLAB 在大学物理热物理学中的一些应用，包括在麦克斯韦速率分布和化工热力学中的应用。

关键词：MATLAB；麦克斯韦速率分布；热力学

Application of MATLAB in thermal physics Abstract:Based on MATLAB’s features of numerical calculation, visualization of graphics processing,opening and scalable architecture,introduced the applications of language of Matlab with high-performance in thermal physics of university physics,include in the Maxwell speed distributionand the chemical industry thermodynamics.

Key Words: MATLAB; Maxwell speed distribution; thermodynamic

引言

热物理学是化工研究的一个方向, 由于热物理学处理比较复杂,恰当地使用可视化以展现数学公式的物理图像, 使其变得直观、形象。MATLAB 是一套高性能的数值计算和可视化软件,下面我们从一些典型的实例出发, 介绍 MATLAB 在热物理学方面的具体应用。

1 MATLAB-PDEtool 介绍

MATLAB-PDEtool 提供了一个功能强大使用灵活的二维有限元偏微分方程求解环境，其图形用户界面更是使用十分方便、直观一般来说，MATLAB-PDEtoo l 包括3个步骤：

（1） 定义一个PDE 问题，它包括确定二维求解区域、边界条件和PDE 系数。MATLAB-PDEtool 能够求解的PDE 型式有：椭圆型、抛物线型、双曲线型、特征值型。当使用GUI 时，可以在画图模式下确定求解区域；在边界模式下选择方程形式和设置方程系数。

（2） 数值求解，它包括剖分、离散方程和得到一个数值解。在GUI 中，在剖分模式下形成满意的网格；在求解模式下通过选择数值计算方法求解。

（3） 图形化显示结果。通常用于表现有限元计算结果的图形有：变形网格图、云图、等值线图、矢量图、网格图、表面图、流线图等。

2 MATLAB 在麦克斯韦速率分布中的应用

气体动力学理论中麦克斯韦速率分布律是大学物理讲授与学习中的一个难点和重点。这是因为公式比较复杂抽象，数学推导证明比较繁琐。如果借助 Matlab 就可以比较方便地解决这些问题。首先, 推导三种速率和归一化条件。已知分布函数表达式为：

KT mV v e v kT

m f

222

3)()2(4-=ππ （1） 最大概然速率分布可由下式求出：

0)(=dv

df v （2）

平均速率的定义是：

dv vf v o

v ?

∞

=)(

（3）

方均根速率为： 21)

(0

2)(2

dv f v v rms

v v

?∞==

（4） 归一化条件是：

dv vf v ?

∞

)( （5）

后三项求解比较复杂,其中用到Gamma 函数,传统方法是查数学用表得到结果。如果应用 Matlab 的符号计算功能, 只需要简单几行语句就可以解决这些问题。用到相关函数有：符号变量创建函数 syms ，求微分函数diff ， 求积分函数 int ，符号化简函数 simple ，字符串转化函数 eval 。

具体 Matlab 语句如下：

Syms mktvanp

g=exp(- (m*v^2)/(2*k*t))*v^2；

a=m/(2*k*t);f=int(exp(- a*v^2)*v^n,v,0,inf)； b=f*4*pi*(a/pi)^(3/2)；

pretty(solve(diff(g,v)))； %求解最概然速率 n=2;eval(simple(eval(b)))，%证明归一化条件 n=3;pretty(simple(eval(b)))，%求解平均速率

n=4;pretty(eval(simple(eval(b^0.5))))， %求解方均根速率 程序中第五行语句给出最大概然速率：

m

k T

v

p

2=

（6） 其物理意义是：“若把整个速率范围分成许多相等的小区间，则 vp 所在的区间的分子数占分子总数的百分比最大”。所以可知在( 0，+∞) 速率区间的分子数占分子总数的百分比(对应着曲线与X 轴所围的面积) 恒等于1，即满足归一化条件，第六行语句证明了这个结论。

第七行语句给出平均速率：

m

kT

v 8=

（7） 第八行语句给出方均根速率：

m

kT

v

rms

3=

（8） 其次,运用 Matlab 强大的画图功能，可以画出同一种气体分子( 氮气) 在不同温度下的分布曲线( 图 1)和不同气体分子( 二氧化碳、氧气、甲烷) 在同一温度下的曲线( 图 2)。程序中用到画图函数 ezplot 。

图1、图2中实线为麦克斯韦速率分布曲线，虚线的 X 轴坐标为最大概然速率 p v 的值。从图中曲线变化情况可以验证p v 是气体分子质量 m 的减函数， 是温度 T 的增函数。由此可以说明在满足归一化条件下，温度升高时曲线变得平坦些，并向高速区域扩展。即温度越高，速率大的分子越多，这就是通常所说的温度越高，分子运动越剧烈的真正含义。

3 Matlab 在处理化工热力学数据中的应

3.1 利用Matlab 处理实验数据的原理及方法

非线性最小二乘法的数学模型为:

)(min 12x f m

i i ∑==

（9） 其中：),...,,(321n x x x x x =是n 维向量，,))(),...(),(()(21T m x f x f x f x f =则上式又可 表示为 22||)(||min x f

（10） 求解非线性最小二乘问题的算法主要有Gauss-Newton 法和Levenberg-Marquadt 法。Gauss-Newton 法通过在每一次迭代步骤求解下列线性最小二乘问题来获得搜索方向k d 。

22||)()(||min k k k x F d x J -

（11） 搜索方向k d 可以用于一维搜索,以保证每次迭代都使f(x)减小。Levenberg-Marquadt 法(又称阻尼牛顿法)则是用下式求搜索方向k d 。

)()())(((k k k k T

k x F x J d I x J x J -=+λ （12）

其中k λ为阻尼因子。它可以控制k d 的大小和方向。当k λ=0时,即为Gauss-Newton 法；当k λ→∞时，即为最速下降法。

在Matlab 的优化工具箱中有求解非线性最小二乘问题的函数Lsqnonlin():[x,resnorm,residual,exitflag,output,lanbda,jacobian]=Lsqnonlin(‘fun’,x0,LB ，UB,options,p1,p2,…)。在这个函数中，输入变量中只有‘fun ’和0x 是不可缺省的，fun 是给出向量函数的m-文件的文件名，0x 是x 的初值。输出宗量中只有x 是不可缺省的，

它是问题的解。针对不同的问题，可以选择不同的算法来计算。默认时，该函数将采用Gauss-Newton 法进行一维搜索，在残差较小时，这种算法有较好的效果和较快的速度。也可通过将输入变量https://www.wendangku.net/doc/763104209.html,rgescale 设置为‘off’来使用Levenberg-Marquardt 算法进行一维搜索[11]。下面我们就来运用这一函数来寻求最佳配偶参数1A 和2A ，并将得到的参数代入方程(1)来检验计算值和

实验值的误差。

(1)建立以汽相组成误差的平方和作为目标函数:

])()[(2222114

11j c e j c j e y y y y F -+-=∑= （13）

其中1e y 和2e y 为实验值，1

c y 和2c y 为待计算值。

(2)建立jscs.m 文件:

%计算配偶参数

Function z=jscs(a)

x1=[0.0000,0.0098,0.0495,0.0811,0.1020,0.1497,0.2494,0.3213,0.4501,0.6905,0.9112,0.9656,0.9818,0.9999];

x2=1-x1;

y1e=[0.0000,0.0148,0.0603,0.0951,0.1102,0.1454,0.1988,0.2271,0.2625,0.3240,0.5551,0.7953,0.8999,0.9999];

y2e=1-y1e;

p=[101.325,101.325,101.325,101.325,101.325,101.325,101.325,101.325,101.325,101.325,101.325,101.325，101.325,101.325];

ps1=[43.2649,42.6417,41.7209,41.4176,41.2667,41.2667,41.5690,42.0260,43.4219,49.7470,71.2944,85.8818,92.1282,101.2401];

ps2=[101.7462,100.1531,97.8025,97.0291,96.6444,96.6444,97.4152,98.5809,102.1477,118.4094,174.7990,213.6398,230.4018,254.9750];

for i=1:14

y1c(i)=exp(a(1)*x2(i)/(x1(i)+a(1)*x2(i)))*ps1(i)*x1(i)/(p(i)*(x1(i)+a(1)*x2(i))*exp((a(2)*x2(i)/(x2(i)+a(2)*x1(i)))));

y2c(i)=exp(a(2)*x1(i)/(x2(i)+a(2)*x1(i)))*ps2(i)*x2(i)/(p(i)*(x2(i)+a(2)*x1(i))*exp((a(1)*x1(i)/(x1(i)+a(1)*x2(i)))));

z(i)=(y1e(i)-y1c(i))^2+(y2e(i)-y2c(i))^2; end

%计算平均偏差

y1c=(ps1.*x1.*exp(a(1).*x2./(x1+a(1).*x2)))./(p.*(x1+a(1).*x2).*e

xp(a(2).*x2./(x2+a(2).*x1)));

n=size(ps);

meandev=sum(1.0-y1c./y1e)/n(2) (3)在Matlab 命令窗口输入:

a =[1 1]；

%0a 为参数1A ,2A 的初值 x=lsqnonlin(‘jscs’, 0a ) 运行后,具体结果见表1。

表1 计算结果

序号 1e y

1

c y

c y （文献计算值）

序号 1

e y 1

c y c y （文献计算值）

1 0.00000 0.00000 0.00000 9 0.26250 0.27690 0.23943

2 0.01480 0.01810 0.00861 10 0.32400 0.37930 0.37318

3 0.06030 0.07580 0.04032 11 0.55510 0.64730 0.67347

4 0.09510 0.10960 0.06253 12 0.79530 0.81960 0.84157

5 0.11020 0.12780 0.07603 13 0.89900 0.89300 0.90940

6 0.14540 0.16170 0.010395 14

0.99990 0.99910

0.19697 7 0.19880 0.21010 0.015338 平均偏差 -0.10093

0.19690 8

0.22710

0.23590

0.18481

配偶参数

1A =0.3994

2A =0.5579

注：平均偏差=[

∑-14

1

0.1(计算值÷实验值]]÷14

4.2 小结

汽相组分的计算值与实验值的误差能够满足实验的要求，达到了计算的目的，说明matlab 的计算程序是有效可行的。本文得到的计算值与文献的计算值相比误差更小，并且编程简单，体现了matlab 进行科学计算的优越性。

4 总结

由以上论述可以看出MATLAB在热物理学方面有很大的用处，方便快捷清晰形象的解决了热物理学中一些复杂的问题。在实际应用中，更多的复杂问题都可以用MATLAB方便快捷的解决，应该根据具体情况，灵活运用MATLAB工具。

5参考文献

1.石玉芳，陈静，刘跃进.MATLAB在化工模拟计算中的应用.[J] 计算机与应

用化学2005,22（9）.

2.李梦龙，王智猛，姜林，刘丽霞.化学软件及其应用. [M]化学工业出版社

2004.10.

3.钱蕙，王锡森. MATLAB在物理化学实验数据处理中的应用. [J] Journal of SuZhou Railway Teachers College 2002,19（1）.

4.李克娥，苏铁军 .MATLAB在处理化工热力学数据中的应用.重庆科技学院学报（自然科学版）.2005,7（3）.

5.李本详，董新荣，曾盔.MATLAB在化工数据及模拟参数拟合中的应用.计算机与应用化学.2007,24（9）.

6.沈春银.MATLAB在化工实验研究中的应用》.南通工学院学报（自然科学版）.2002,1（1）.