微分计算

2024-09-27

微分计算（共8篇）

微分计算篇1

1计算函数增量的近似值

例1半径10, 厘米的金属圆片加热后半径伸长了0.05, 厘米问面积增大了多少?

解:设A=πr2, r=10厘米, ∆r=0.05厘米.

∴∆A≈dA=2πr⋅∆r=2π×10×0.05 2=3.14 (厘米) .

例2:设有边长为1m的立方体, 其表面要涂上一层厚为0.01cm的贵重金属, 该贵重金属的比重为10g/cm3。试估计一下需多少克金属?

解:设立方体的体积为V, 要计算当边长l由1m=100cm增加到100.01cm时, 体积V增量是多少, 因为∆l=0.01cm很小, 所以可以用d V代替∆V, 而计算d V比计算∆V简捷些.为此要找出V与l的关系, 易知

2计算函数的近似值

3误差估计

由于测量仪器的精度, 测量的条件和测量的方法等各种因素的影响, 测得的数据往往带有误差, 而根据带有误差的数据计算所得的结果也会有误差, 我们把它叫做间接测量误差。

定义:如果某个量的精确值为A, 它的近似值为a, 则叫做a的绝对误差。

而绝对误差与的比值叫做的相对误差。

问题:在实际工作中, 绝对误差与相对误差如何求得?

方法:将误差确定在某一个范围内。

如果某个量的精确值A是, 测得它的近似值为a, 又知道它的误差不超过, 即

则, 叫做测量的绝对对差限而叫做测量A的相对误差限.

通常把绝对误差限与相对误差限简称为绝对误差与相对误差.

例5设测得圆钢截面的直径D=60.03mm, 测量D的绝对误差限利用公式计算圆钢的截面面积时, 试估计面积的误差.

本文主要通过微分的概念, 总结出微分的三种应用, 为解决实际问题提供有力的帮助。

摘要：微分在近似计算中的应用主要有三方面的内容:计算函数增量的近似值;计算函数值的近似值;误差估计。

关键词：微分,近似计算

参考文献

[1]数学分析 (第二版) [M].华东师范大学数学系.

[2]盛祥耀.高等数学辅导[M].

微分计算篇2

关键词：压缩映射;不动点;微分;泰勒公式;近似计算

Abstract：We are difficult to calculate the accurate numbers of the results in production and life and scientific research，but can use approximate numbers to replace them within the limits of allowable error. In this paper，three methods of application for compress mapping theorem and differential and Taylor formula in approximation algorithm are showed.

Keywords：compress mapping;fixed point;differential;Taylor formula;approximation algorithm

一、压缩映射原理在近似计算中的应用

1 压缩映射原理

定义（泛函分析中的定义）给定度量空间及映射，若存在，使得，则称为映射的不动点，若存在正常数，使得对任给的，都有成立，则称是上一个压缩映射.

定理（巴拿赫压缩映射原理）设为完备的度量空间，是压缩映射，则有且仅有一个不动点.

注压缩映射是连续的.

注空间的完备条件是为了保证映射的不动点存在.

定义（高等数学上的定义）设在闭区间上有定义，方程在上的解被称为在上的不动点.若存在常数，使得，都有成立，则称是上的压缩映射.

定理设是有界闭区间上的压缩映射，那么函数在闭区间上存在唯一的不动点.

推论如果函数在闭区间上可导，且，则函数在闭区间上存在唯一的不动点.

推论如果函数是到自身的映射，在上连续，且，则函数在上存在唯一的不动点.

定理设函数是闭区间上的压缩映射，且，而，，对，有，则函数在上存在唯一的不动点，而且 .

注若非闭，函数是压缩映射，但是函数在中的不动点未必存在.

2 在近似计算中的应用

例1 试用压缩映射的原理计算的近似值.

解注意到是方程

的是实根，构造辅助函数

，

则任给，都有

，

可令

，

容易验证，当时，有

，

所以是压缩因子的压缩映射.由于是中的有界闭集，故有，使得进而可用迭代法求得的近似值.

取，从而有

，

，，，

由上述说明可知

，

从而还可求出近似值与精确值之间的误差.

由上可知，巴拿赫不动点定理不仅能够证明一定条件下不动点的存在性和唯一性，而且提供了计算不动点的方法，即迭代法.一般地，我们可以从任意选取的一点初值出发，逐次作点列的迭代运算，其最终收敛到所求方程的解.这种方法又称为逐次逼近法，这也是计算数学中的一种重要方法.

二微分公式在近似计算中的应用

1 微分的概念

定义设函数在某区间内有定义，当自变量在点处产生一个改变量（其中）时，函数的改变量与有下列关系

，

其中是与无关的常数，则称函数在点处可微，称为函数在点处的微分，记为 .

注4 由可微与可导的关系，进一步可得

忽略掉高阶无穷小，有 .

2 在近似计算中的应用

例 2 试用微分计算的近似值.

解注意到，且可直接开方得2，于是设函数

，

有

，且，

取给一个增量，对应函数增量

即 .

例 3 试用微分计算的近似值

解注意到，且，于是设函数

，

有

，

取给一个增量，对应函数增量

即 .

由于近似公式里，省略了高阶无穷小，因此，选择微分作近似计算时，误差取决于自变量的增量，随的减小而减小.

三泰勒公式在近似计算中的应用

1 泰勒中值定理

定理（泰勒中值定理）如果函数在含有点的某个开区间内具有直到阶的导数，那么对于，有

称为在点处关于的阶泰勒公式.其中（在与之间）稱为拉格朗日型余项.

注5 当时，是比的高阶无穷小，故

注6 在泰勒中值定理中，若取是，泰勒公式变成如下形式

，

称其为的阶麦克劳林公式.

2 在近似计算中的应用

例4 应用阶泰勒公式计算的近似值.

解注意到要计算的近似值，需利用函数在点的解泰勒公式而得到，由阶泰勒公式

得到

误差（在与之间）.

例5 应用阶泰勒公式计算的近似值.

解注意到与0很接近，因此应用函数的3阶麦克劳林公式

得到

误差（在与之间）

由泰勒公式可以观察到，利用泰勒公式作近似计算时，选取泰勒公式的阶数越高，近似计算的精确度越高.并且微分的近似计算公式就是一阶泰勒公式的特殊情况，因此微分近似计算的精确度并不高，但是在精确度要求不高的情况下，选择微分作近似计算更加简单.

参考文献：

[1] 刘炳初. 泛函分析[M]. 北京：科学出版社. 1998：3-22.

[2] 邓志颖，潘建辉. 巴拿赫不动点定理及其应用[J]. 高等数学研究，2013（4）：78-81

[3] 强文久，李元章，黄雯荣.数学分析基本概念与方法[M]. 北京：高等教育出版社，1993：29-35.

[4] 刘玉琏，傅沛仁. 数学分析讲义：上册[M]. 4版. 北京：高等教育出版社. 2003：57-60.

[5] 孔繁亮，张立斌，巨小维，朱捷，高等数学：上册[M]. 2版. 哈尔滨：哈尔滨工业大学出版社. 2011：55-60.

[6] 同济大学数学系. 高等数学：上册[M]. 6版. 北京：高等教育出版社. 2008：139-144.

作者简介：

巨小维（1982-），女，山东临朐人，讲师，主要研究方向：非线性泛函分析。

微分在近似计算中的应用篇3

若函数y=f (x) 在x0处可导, 且f′ (x0) ≠0.当|Δx|很小时, 我们有

上式表明, 若函数f (x) 在x点处的函数值难于计算, 但在x邻近的一点x0处的值容易计算, 且f′ (x0) 的值也容易计算, 则我们可用f (x0) +f′ (x0) Δx作为f (x) 的近似值.如公式 (2) .

另外, 有时要计算函数f (x) 的增量Δy, 由于Δy=f (x0+Δx) -f (x0) 难于计算, 我们也可用f′ (x0) ·Δx作为增量Δy的近似值.如公式 (1) .总之, 公式 (1) 与公式 (2) 是微分在近似计算应用中常用的两个基本公式.

2.应用举例

例1:有一批半径为1厘米的求, 为了提高球面的光洁度, 要镀上一层铜, 厚度定为0.01厘米.估计一下每只球需用铜多少克 (铜的密度是8.9克/厘米3)

解:先求出镀层的体积, 再乘上密度就得到每只球需用铜的质量.

因为镀层的体积等于两个球体体积之差, 所以它就是球体体积V=4/3πR3当R自R0取得增量ΔV.我们求V对R的导数

将R0=1, ΔR=0.01代入上式, 得

于是镀每只球需用的铜约为

0.13×8.9=1.16 (克)

例2:求 (1.0125) 2016的近似值.

解:1.0125自乘2016次难于计算, 我们构造函数f (x) =x2016, 1.0125邻近的一点x0=1, 且f (x0) =f (1) =12016=1.便于计算.并且f′ (x) =2016·x2015,

于是, f′ (x0) =f′ (1) =2016·12015=2016也容易计算, 且此时

从而, 由公式 (2) , 有

例3:计算sin30°30′的近似值

解:把sin30°30′化为弧度, 得

由于所求的是正弦函数的值, 故设f (x) =sinx.此时f′ (x) =cosx.

如果取x0=π/6, 则f (π/6) =sinπ/6=1/2与都容易计算, 并且比较小.应用公式 (2) 便得

3.问题的总结

通过上述三个例子, 我们看到利用微分往往可以较容易地解决一些近似计算的问题.当然, 若在公式 (2) 中令x0=0, 有

利用上式, 我们还可以推出工程上常用的近似公式如下:

(假设|x|是较小的数值) .

参考文献

[1]同济大学数学系.高等数学 (第四版) [M].北京:高等教育出版社, 1996.

微分计算篇4

求解常微分方程组的数值方法不止龙格-库塔法一种, 还有亚当姆斯法等。龙格-库塔法虽然工作量较大, 但是它在计算值时, 只用到, 而不直接依赖于和等。也就是说, 在初值确定后, 就可以依次计算和等等。直至整个过程结束, 不存在计算起步的问题。另外, 这种方法没有规定后一步的步长与前一步步长必须满足的某种关系, 可以任意改变计算步长。龙格-库塔方法每一步须要4次计算, 看起来计算过程烦琐, 实际上在M A T L A B运用程序进行计算非常快速准确。

1 微分方程组的程序框图

编程思路的逻辑框图为图1和图2。

2 微分方程组的仿真计算

由某实例受力及运动分析, 建立起如下方程组:

其中,

把μ1 (Z) 分段函数和以上微分方程组联立, 利用四阶五级龙格-库塔法, 借助MATLAB软件编制微分方程组的求解程序, 对该数学模型进行仿真计算, 得出p, v, l, z关于t的图像, 计算结果如图3、图4和图5。

3 结语

仿真计算结果与实测的p-t、v-t、Z-t曲线一致性较好, 表明该微分方程组的计算正确。利用计算机仿真技术对某数学模型进行计算, 节约了研制经费, 为优化设计提供参考。

摘要：传统的解微分方程组的方法有近似分析解法﹑表解法和图解法。这些方法都有局限性, 电子计算机编码的出现及其应用, 不仅代替了繁重的人工求解, 而且改变了传统的研究方法。MATLAB是一种基于矩阵的数学软件包, 该软件包包括了一个数值程序扩展库, 并且有高级编程格式。应用四阶五级龙格库塔法编制Matlab程序对微分方程组进行求解, 结果表明无论是曲线或是特殊点与试验实测值一致性都比较好

关键词：微分方程组,数值计算,仿真,MATLAB

参考文献

[1]楼顺天, 姚若玉, 沈俊霞.MATLAB 7.X程序设计语言[M].西安电子科技大学出版社, 2007.

微分计算篇5

原子与分子碰撞实验及其理论研究是原子分子物理十分重要的研究方向之一[1], 它为冲击波、声波、风洞流扩张的快速压缩过程中出现的弛豫现象、气相反应和输运性质、气体激光、转动激发的共振荧光过程等科技领域提供了适用的方法和大量的参考信息, 从大量文献[2,3,4,5,6,7,8,9,10]可以看出。

现在计算量大但能更好地处理电子相关能的耦合簇理论CCSD被越来越多地用于研究范德瓦耳斯分子间弱相互作用[11], 如Gerrit[12]等用RCCSD (T) 方法研究了He-O2的三维势能面;李绛[11]等用CCSD (T) 方法研究了Ne-HCl的势能面, 所得结果与实验符合较好。

本文正是基于以上理由, 采用CCSD (T) 方法, 利用大基组和键函数计算得到了He-Al H体系的微分截面, 希望通过大量的计算对该体系的实验研究具有一定的指导作用。

2 计算方法

式中, d'是双原子分子中两核之间的相对位置矢量, o'是入射原子A相对靶分子BC质心的相对位置矢量。体系的Hamilton算符为:

其中

3 计算结果与讨论

对He-Al H体系采用量子化学从头计算的单双迭代, 包含非迭代三重激发微扰的耦合簇CCSD (T) 方法和相关一致基组aug-cc-p V5Z[14,15], 并采用了3s3p2d1f1g基集的高斯键函数。计算He-Al H体系结构模型如图1所示。所有的从头计算均在Gaussian 03程序包中完成, 并进行基组重叠误差 (BSSE) 校正[16], 得到He-Al H体系的相互作用势数据。然后将相互作用能数据构造为以Legendre函数展开的各向异性势解析表达式为

其中, Pλ (COSθ) 是Legendre函数, Vλ (R) 是与R有关的径向系数。He-Al H体系相互作用势能面如图2所示。

根据计算结果, 进一步分析讨论了He-Al H碰撞体系微分截面的散射特征。

图3和图4是仅列出He原子入射能量分别为10, 50, 100me V时, 总微分截面随散射角Θ的变化。从图3和图4中可以看出, 在一定的入射能量时, 总微分截面在散射角Θ=0°时最大;随着散射角增加, 总微分截面呈有规律的散射振荡而逐渐减小, 随后散射振荡消失而逐渐减小。明显的散射振荡区出现在小角区 (Θ<90°) 。随着入射能量的增加, 在散射角Θ=0°的总微分截面值增加;散射振荡的散射角逐渐减小, 散射振幅也逐渐减小。

图5是仅列出He原子入射能量分别为10, 50和100me V时, 非弹性微分截面随散射角Θ的变化。从图5中可以看出, 非弹性微分截面在小角区 (Θ<90°) 也有散射振荡现象, 并随散射角Θ的增加而减小, 散射振幅较小。在大角区 (Θ>90°) , 非弹性微分截面散射振荡现象消失。在较低入射能量 (E=10me V) 时, 非弹性微分截面在大角区 (Θ>90°) 随散射角Θ的增加而增加。在较高入射能量 (E=50, 100me V) 时, 非弹性微分截面在大角区 (Θ>90°) 随散射角Θ的增加而缓慢减小。但在约Θ=20°~50°区间, 非弹性微分截面随入射能量的增加而增加。

4 结论

微分计算篇6

抽水试验是探明含水层水力特征的一个重要手段。通常,抽水试验是采用降深-时间曲线来分析含水层的特征并计算水文地质参数。该种方法的理论依据是地下水向井运动及泰斯(Theis)假设。然而抽水试验要达到稳定要求,需经过较长的持续时间,在这期间降深-时间曲线变化既受到含水系统的内边界,外边界的影响,会受到含水层的非均值性的影响,因此,并不是整个抽水过程都会符合理论计算假设[1]。基于以上的原因,分析判断出抽水试验降深-时间曲线中符合理论假设的有效性的区域,是计算水文地质参数,选择合适模型算法的一个关键性步骤。

基于降深-时间的微分分析方法的抽水试验数据分析,是一个更具代表性的分析方法。该方法较传统的降深-时间曲线对降深的改变有更高的敏感性。这种较高的敏感性有利于用该方法定义出含水系统的内边界条件,确定外边界条件对降深的影响,确定在抽水试验过程中的水流状态及符合地下水向井运动的时间区域。

微分分析法多用于石油工业,用来分析研究油气的运移规律[5]。Tiab和Kumar[6]第一次对定流量抽水实验数据采用了微分分析。随后该法被大量用于地下水流动规律的研究,提高了人们对水力试验数据的认识,明确了不同水流特征的定义。Karasaki,Spane和Sapan及Wurstner[2~4]等人首次将微分分析法引入到对地下水流流动的研究中,对含水层的水流特征进行了分析。Karasaki及Ostrowski和Kloska等人[2,7]在承压含水层中的水位恢复实验使用了微分分析法对数据进行处理。

利用抽水试验数据计算水文地质参数的井流模型及求解方法有很多,主要根据含水层特性分类。如关于承压含水层的数值解法有:Theis(1935)/Hantush(1961)解法,Theis(1935)解法,Cooper-Jacob(1946)解法,Moench-Prickett(1972)解法等,关于潜水含水层的数值解法包括:Neuman(1974)解法,Moench(1997)解法,Artakovsky-Neuman(2007)解法等,关于裂隙含水层包括:Moench(1984)解法,GringartenWitherspoon(1972)解法,Gringarten-RameyRaghavan(1974)解法,Gringarten-Ramey(1974)解法[8]。但大多未对降深-时间数据进行处理,直接配线利用公式计算参数。对数据采用微分分析的,并根据微分分析结果选择合理求解的方法,并不多见。

1 抽水试验数据的微分分析

1.1 研究区域水文地质条件

研究区域内主要的含水层为石炭系中上统壶天群(C2+3ht)岩溶含水层,其岩性为白云岩、白云质灰岩,广泛分布于研究区域内。受褶皱和断层错动的影响,含水层底板标高变化非常大,其平均厚度约为149.37m,平均顶板标高为90.04m,天然状态下静止水位标高平均为101.36m,具有微承压性。

壶天群(C2+3ht)含水层上覆一层10~30 m的第四系粉质粘土、粉砂土,由于岩溶发育,岩溶含水层与第四系含水层有着紧密水力联系。研究区北部及东北部有石炭系下统测水组(C1)和泥盆系上统帽子峰组(D3m)砂页岩隔水层出露,研究区西部有泥盆系上统天子岭组中、上亚组(D3tbc)杂质灰岩相对隔水层分布,组成了北、西隔水边界。研究区南部及东南部壶天群白云岩分布广泛,具备本含水层的地下水向研究区径流补给的边界条件。图1为研究区域的水文地质简图。

壶天群(C2+3ht)含水层岩溶发育,钻孔岩溶率为5%左右,溶洞大小一般为0.5~3m,其长度较短,不超过几米,洞与洞之间靠裂隙联接,且溶洞多为全填充溶洞。含水层富水性较强,不同地段的钻孔单位涌水量从0.03~0.06L/s.m。对于同一水文地质单元的壶天群含水层,尽管岩溶发育存在着明显的不均一性,但在强岩溶带内水力联系较密切,宏观上有统一的地下水水位。

1.2 微分分析模型

降深-时间数据的分析,是解译抽水过程含水层特征的重要部分。抽水试验进行过程中,水头会随时间不断降低,通过降深-时间曲线可以计算出含水层的特性。经典的Theis解法,采用双对数坐标,该方法使得时间和降深在坐标中的分布不一致。根据Cooper-Jacob[9]解法绘制的半对数降深-时间曲线,尽管克服了降深-时间分布不一致的特点,但在该曲线中,很难区别出哪段曲线是满足于Theis假设及无限含水层向井运动假设的。由于微分分析中,时间对降深有较强的敏感性,使抽水过程受边界条件及含水层非均质性影响的部分在降深-时间曲线中有了较明显的区分。图2,对比了采用双对数的降深-时间曲线与采用微分分析处理后的降深-时间曲线,列出了理想的降深微分分析曲线的不同水流形态和边界条件。

降深的微分分析,考虑降深随着时间对数的变化,由于降深的变化是一个瞬时的变化,不能够直接测量得到,因此考虑一个平均值。降深导数斜率的计算,使用以下的公式[10]。

其中,ti,ti+j和ti-k分别为,斜率的中心所对应的时间点,比ti早0.1~0.5时间点,比ti晚0.1~0.5个时间点;si,si+j和si-k,分别为ti,ti+j及ti-k时间点所对应的水位降深。

1.3 抽水试验微分分析

抽水试验微分分析数据采用壶天群(C2+3ht)含水层,历时24h的定流量抽水试验,抽水流量为3.325L/s,最大降深达11m。停泵后,恢复水位历时200min。

图3~图5描述了抽水井的抽水试验结果,并对比了采用半对数坐标的降深-时间的曲线(图3),双对数坐标的降深-时间曲线(图4)及微分处理后的降深-时间曲线(图5)的结果。从图中可以看出,图5所表达的抽水过程,降深随时间的变化,存在明显的三个阶段。而图3与图4,降深随时间虽有变化,但不明显,三个变化阶段不易区分。

分析图5中三个阶段,分别代表了重力释水的延迟,无限向井运动,边界对抽水试验的影响。2~190 min这段时期,降深随时间变化形成一个峰谷,这是一个过渡时期,受重力释水的延迟作用,通常会出现在抽水消耗完井中原有储量及形成无限向井水流之间。在抽水试验190~400min这段时期内,存在一个斜率接近0的直线段,该段降深-时间曲线表现出了无限向井运动的水流特征。400min后直至实验结束,该段曲线几乎垂直于坐标横轴,降深随时间不发生变化,可以认为该段是由于存在一个定水头边界所致。由于0~2min之内的降深数据未能及时收集到,该段降深变化趋势未有明显的规律,并不能判断井中已有水量是否对抽水试验的前期数据造成影响。

2 水文地质参数计算

2.1 井流模型确定

由以上对抽水试验时间-降深曲线的微分分析(图3~图5),表明该抽水过程,有两个较为明显的阶段,水位降深的延迟期及径向流期。这一特征与潜水含水层在抽水阶段所具有的典型的特征较为一致,采用考虑流速垂直分量和弹性释水的纽曼模型(公式(2)~(7)),

其数学描述为:

其中,Kr为水平径向渗透系数;Kz为垂向渗透系数;μs为储水率;μ为给水度;H0为潜水流初始厚度。

2.2 水文地质参数计算

利用配线的方法,根据公式(2)~(7),可以得到含水层的水文地质参数。其计算结果如图6及表1所示。

从图6中可以看出由微分分析判断选取的含水层类型理论的降深-时间曲线与实测的降深-时间曲线匹配度较高,在参数计算时选取试验后期数据是位于径向流区域的数据,其计算出的水文地质参数也更加可信。

根据计算的参数(表1),可看出试验区域内含水层较一般的空隙含水层具有较高的给水度及导水系数,且渗透性在垂向及纵向上有着明显的差异,其纵向渗透系数与横向渗透系数(Kv/Kh)比值为4.01,因此可认为垂向渗透性对含水层的水文地质特征有着较大的影响,这一特征反映出一般浅层岩溶含水层岩溶纵向发育特点,这不同于孔隙介质;降深-时间曲线中(图6),出现一段水位降深的滞后期,由于所测水头,均高于地面,有一定的承压性,是造成降深滞后的原因,主要是受到上覆第四系含水层影响,反映出他们之间有着较强的水力联系。图6水位降深的延迟期与理论曲线还存在一定的偏差,试验降深略大于理论降深。是因为试验含水层为岩溶含水层,尽管它具有统一的水力面,但由于长期受岩溶作用得影响,含水层的非均质性表现强烈,会造成这一偏差,因此在利用孔隙介质的理论假设去模拟岩溶介质时应该要更为慎重。

3 结论

(1)由于岩溶含水层表现出较强的非均质性特性,因此在抽水试验过程中区分出受不同条件影响的降深变化区段,选取符合理论假设的区段进行参数计算是分析抽水试验数据判断含水层特性的一个重要步骤。

(2)由于微分分析方法对水头的降深有着较高的敏感性,在整个抽水过程中,水头降深受到内外边界影响而表现出的差别,在微分分析中更加的明显,因此利用该方法来分析岩溶含水层的水文地质特征更加客观和有效。

微分计算篇7

现代复杂机电产品（如航空航天器、机器人、汽车等）通常是集机、电、液、控、磁等多学科领域于一体的复杂物理系统，经常表现出时间依赖（对时间导数）与空间依赖（对坐标偏导）共存的行为特征，而且可能呈现出多领域之间及时间域与空间域之间的耦合特性[1]。

物理系统行为规律的描述通常有两种主要的形式。系统在单纯时间域的物理行为往往由常微分方程（ordinary differential equation,ODE）描述，如果涉及代数约束，则形成微分代数方程（differential-algebraic equation,DAE),DAE是描述时间域物理规律的普遍形式，如机械多体、电子电路等系统规律的描述；若物理行为涉及空间场，出现对空间变量的偏导，则往往由偏微分方程（partial differential equation,PDE）描述，如位势、传热、波动等相关物理系统的规律描述[1]。

物理系统建模经历了从面向过程建模到面向对象建模、连续域与离散域分散建模到连续-离散混合建模、单一领域独立建模到多领域统一建模的发展阶段。Modelica语言是近年来欧洲仿真界为解决复杂物理系统建模与仿真问题而提出的一种多领域统一建模语言[2,3]，然而，目前Modelica语言只能对时间域的物理系统进行统一建模，还不能对空间域的物理系统进行描述，更无法对其进行仿真优化，这大大限制了Modelica语言的应用范围。为此，国外已有学者着手扩展Modelica语言以支持PDE问题的建模与仿真[4]。周凡利等[1]提出了解决该问题的思路，但没有具体实现。李志华等[5,6]也开展了这方面的研究工作，初步实现了Modelica语言对PDE问题的描述、建模以及仿真求解。然而现有的这些方法都是采用简单的有限差分法或线上法来求解具有规则区域（如矩形）的PDE问题，对于不规则区域的复杂PDE问题则无法求解。

本文在李志华等原有工作[5,6]的基础上，从多领域统一建模与仿真的角度，针对一般性的PDE问题（包括不规则求解区域、复杂边界条件、线性或非线性PDE系统），提出一种PDE与DAE的一致求解方法，为Modelica直接求解复杂工程系统中多领域耦合、时间域与空间域耦合的复杂问题奠定基础。

1 PDE的已有解法

PDE的解法主要分为解析法和数值法。到目前为止，只有有限形式的PDE能够得到解析解，在工程实际中一般采用数值求解。PDE的数值求解技术比较成熟，可用的方法包括有限差分法、有限元法、有限体积法以及线上法等[7]。但有限差分法和线上法只适用于规则的求解区域，而有限元法和有限体积法是基于网格的计算方法，在某些工程问题（如动态裂纹扩展、高速撞击、冲击破坏、流固耦合等）中存在网格的束缚，使得计算遇到很大的困难，因此出现了无网格方法[8]。

无网格方法只需要节点的信息，不需要节点与节点之间相互联系的信息，这样很容易在复杂计算区域内布置节点。无网格方法的构建主要包括近似函数的构建方法和微分方程的离散方法两个部分，根据近似函数构建方法和微分方程离散方法的不同，可以构建出许多不同的无网格方法[8]。目前比较常用的近似函数的构建方法有：核函数近似方法、再生核近似方法、移动最小二乘近似方法、单位分解近似方法、径向基函数近似方法、点插值近似方法等。微分方程的离散方法包括加权残量法、配点法、Galerkin法以及局部Petrov-Galerkin法。

Khattak等[9]运用无网格配点法成功求解了一类非线性PDE;Yao等[10]应用全局和局部径向基函数来求解三维抛物线型PDE，并比较了这两种无网格方法的性能；Kamruzzaman等[11]利用多项式点插值和配点法来构造无网格方法，较好地求解了椭圆形、抛物线型和双曲线型PDE；吴宗敏[12]介绍了散乱数据拟合研究中的径向基函数方法，及其在散乱线性泛函信息插值、无网格PDE数值解中的应用；吴孝钿[13]应用Sobolev splines径向基函数和紧支柱正定径向基函数，得到求解PDE边值问题的无网格算法，并针对散乱数据的特点，给出了计算整体稠密度h的算法及如何通过加密节点使h值缩小的一个可行方法。然而上述无网格方法都是直接对时间变量和空间变量进行离散，只适合求解单纯的PDE问题，不适合求解复杂的PDE与DAE耦合问题。

本文借鉴传统的求解PDE的无网格方法，选用径向基函数和配点法来构建径向基函数配点无网格法，并对其进行改进，即只对空间变量离散而保持时间变量连续，直接将PDE在空间上（即配点处）离散成一系列的DAE，然后利用成熟的DAE求解器进行统一求解。

2 径向基函数配点无网格法

对于d维实空间中定义在域Ω上的PDE问题：

式中，L、B为微分算子；u(X,t）为未知量场函数；f(X,t）、g(X,t）为已知函数。

我们所构建的径向基函数配点无网格法的基本思想是：首先采用径向基函数构造近似函数，将未知量场函数的时-空变量分开，然后运用配点法对空间变量进行离散，而保持时间变量连续，这样就将PDE问题在空间上离散成一系列只含时间变量的DAE问题，具体过程如下。

首先在PDE的不规则求解区域Ω内和边界Ω上选定N=Nu+Nb个离散的节点（即配点），然后应用径向基函数构造u(X,t）的近似函数，并将其构成时间与空间分离的形式：

其中，N为节点总数；Nu为域内节点数；Nb为边界节点数；αi(t）为待定系数；Xi为实空间上的配点；φi（‖X-Xi‖）为径向基函数，φi（‖X-Xi‖）=φ（ri(X));ri(X）为Euclidian范数，ri(X）=‖X-Xi‖。

为确保解的唯一性，（ri(X））必须无条件正定，这种径向基函数包括Gaussian函数e-cr2、逆MQ函数（r2+c2）β（β<0）和紧支正定径向基函数等。本文采用Gaussian函数。

将式（2）代入式（1）中，并使这N个点满足式（1）的微分方程和边界条件，得到

当微分算子L、B为空间变量的偏导时，由于空间变量已与时间变量分离，且径向基函数φi（‖X-Xi‖）对空间变量可求导，在每一节点处，空间坐标已知，因此Lφ（ri(Xk））或Bφ（ri(Xk））就是一个已知值，这样，式（3）中就只剩下时间的导数，即每个节点处对应一个只与时间有关的DAE，这样就将PDE转化为一系列的DAE（具体过程可参见实例部分）。

进一步，PDE问题（式（1））可变为求解一个N×N的线性方程组问题，用矩阵表示为

其中，S为N×N的矩阵，S=(Ski)N×N;a为待求的系数向量，a=（α1(t），α2(t），…，αN(t));b=(b1,b2，…，bN）；且

由式（4）求出未知系数αi(t)(i=1,2，…，N）后，通过式（2）就可以获得域内和边界上任意一点的场函数值u(X,t）。

3 实例及求解过程

通过编程，利用径向基函数配点无网格法对PDE进行空间离散，自动将PDE问题转化成DAE问题。本文采用MATLAB编程，首先将带时间域的PDE问题进行时间与空间分离（若不带时间域则不用分离），然后通过径向基函数配点无网格法对空间变量进行离散，得到一系列离散点处的DAE，并把这些DAE数据用mat格式保存起来，接着将该mat格式文件导入到基于Modelica语言的多领域统一建模与仿真平台MWorks中[14]，并利用其成熟的DAE/ODE求解器进行PDE与DAE的一致求解。整个流程如图1所示。

以下面一个不规则区域的二维热传导问题为例来说明本文所提方法的有效性：

式中，T为某个时间值。

其求解区域由如下边界组成：

对该不规则求解区域用配点法进行不规则离散，得其节点分布如图2所示，其中域内节点数Nu=61，边界上节点数Nb=42。然后对这些节点从左到右、从下到上进行编号。

对该二维热传导PDE问题，运用上述PDE与DAE的一致求解方法和过程进行一致求解，令

对于求解域内的Nu个离散节点（xi,yi)(i=1,2，…，Nu），满足域内的偏微分方程，即

对于求解域边界上的Nb个离散节点（xi,yi)(i=1,2，…，Nb），满足边界条件方程，即

对于域内及边界上的所有离散节点（xi,yi(i=1,2，…，N),N=Nu+Nb，满足初始条件方程，即

本文采用Gaussian径向基函数，在各离散节点处均可计算出它们的值。这样，式（5）～式（7）中就只剩下时间的变量，即利用径向基函数配点无网格法已将PDE在离散节点处转化为一系列的DAE，然后就可以在MWorks环境中利用其自带的DAE求解器进行求解，从而方便地得出场函数在各个离散节点处随时间变化的函数值。图3表示的是场函数在编号为15节点处的仿真结果；图4所示为本文的数值解与其精确解u(x,y,t)=e-2tsin(x+y）之间的比较，此处t=0.5s。

由图4可以看出，本文的数值解uP与精确解uQ非常吻合，达到了较高的求解精度。进一步，定义本文的数值解与精确解之间的误差为，通过计算得到er=7.1×10-5。

4 求解精度影响因素分析

为了更好地应用径向基函数配点无网格法来求解PDE问题，本文研究了不同离散节点数、平均节点间距和径向基函数参数c的取值对求解精度的影响，如表1所示。

由表1可以看出，一般情况下，当参数不变时，离散节点数越大、平均节点间距越小，则求解精度越高。例如，在c=1不变的情况下，当节点数为14、平均节点间距为0.47时，误差为2.2475×10-4；而当节点数为30、平均节点间距为0.28时，误差则为1.7142×10-5。同时还可以看出，随着参数c的不断增大，求解精度并不是呈递增或递减状态，而是有起伏变化，只有当c取适当的值时，误差er才较小。由此可见，恰当确定径向基函数参数c的取值很关键，然而目前还没有规律可循，只能通过多次反复运算来确定一个合适的值。

5 结论

多领域统一建模要求用偏微分方程和微分代数方程来统一描述和统一求解，本文针对一般性的偏微分方程问题，提出了偏微分方程与微分代数方程的一致求解方法，给出了该方法的实现过程，分析了离散节点数和径向基函数参数对求解精度的影响，得到如下结论：

(1）与传统的无网格方法相比，本文采用只对空间变量离散而保持时间变量连续的策略，能方便地将偏微分方程在离散节点处转化为一系列微分代数方程，从而在不改变Modelica语法的前提下，较好地实现偏微分方程与微分代数方程的一致求解，大大简化了复杂的偏微分方程与微分代数方程耦合问题的求解难度。

(2）实例结果表明，本文所提出的方法能较好地解决具有不规则求解区域的偏微分方程问题，且求解精度高，这有利于Modelica直接求解复杂工程系统中多领域耦合、时间域与空间域耦合的复杂问题。

摘要：Modelica语言是一种复杂物理系统多领域统一建模语言,但目前该语言只能解决由微分代数方程(DAE)描述的问题,而不能解决由偏微分方程(PDE)表达的问题。为此,提出一种偏微分方程与微分代数方程的一致求解方法,利用所构建的径向基函数配点无网格法直接将偏微分方程在空间上离散成一系列的微分代数方程,然后采用成熟的微分代数方程求解器进行求解。实例结果表明,该方法在不改变Modelica语法的前提下,能较好地实现偏微分方程与微分代数方程的一致求解,且求解精度高、边界条件处理简单,有利于Modelica直接求解复杂工程系统中多领域耦合、时间域与空间域耦合的复杂问题。

微分方程应用举例篇8

1 应用问题举例

1.1 生态系统中的弱肉强食问题

在这里考虑两个种群的系统, 一种以另一种为食, 比如鲨鱼 (捕食者) 与食用鱼 (被捕食者) , 这种系统称为“被食者—捕食者”系统。

Volterra提出:记食用鱼数量为x (t) , 鲨鱼数量为y (t) , 因为大海的资源很丰富, 可以认为如果y (t) =0, 则x (t) 将以自然生长率r (r>0) 增长, 即.x=rx。但是鲨鱼以食用鱼为食, 致使食用鱼的增长率降低, 设降低程度与鲨鱼数量y (t) 成正比, 于是相对增长率为r-λy。常数λ>0, 反映了鲨鱼掠取食用鱼的能力。如果没有食用鱼, 鲨鱼无.法生存, 设鲨鱼的自然死亡率为d, 则y=-dy。食用鱼为鲨鱼提供了食物, 致使鲨鱼死亡率降低, 即食用鱼为鲨鱼提供了增长的条件。设增长率与食用鱼的数量x (t) 成正比, 于是鲨鱼的相对增长率为-d+µx。常数µ>0, 反映了食用鱼对鲨鱼的供养能力。所以最终建立的模型为:

这就是一个非线性的微分方程。

1.2 雪球融化问题

有一个雪球, 假设它是一个半径为r的球体, 融化时体积V的变化率与雪球的表面积成正比, 比例常数为k>0, 则可建立如下模型:

其中, 带入上式得到如下一阶线性微分方程:

1.3 冷却 (加热) 问题

牛顿冷却定律具体表述是, 物体的温度随时间的变化率跟环境的的温差成正比。记T为物体的温度, Tm为周围环境的温度, 则物体温度随时间的变化率为, 牛顿冷却定律可表示为:

其中k是正的比例系数, 牛顿定律表示冷却过程时,

且:

T>Tm;

表示加热过程时,

且:

看到, 牛顿冷却定律也是一个一阶线性微分方程。