非线性计算不稳定

2025-01-31

非线性计算不稳定（精选3篇）

非线性计算不稳定篇1

一、引言

根据物理定律, 比如牛顿第二定律、质量守恒定律、能量守恒定律、气体试验定律等, 可以得到支配大气运动的基本方程组。但由于大气运动基本方程组是一组高度非线性的偏微分方程组, 很难求得其解析解, 人们可通过数值的方法求其近似解 (数值解) , 这就是数值天气预报。数值天气预报是一门实用性很强的应用基础学科[1]。通过差分方法求解大气运动基本方程组时, 人们发现数值积分过程中会产生计算不稳定问题, 这就需要采用恰当的差分格式或积分格式。在《数值天气预报》课程关于时间积分格式的讲解中[2], 重点介绍了“显示”格式, 即差分方程的右端项全部为当前 (或和过去) 时刻的变量值, 通过积数值分可求出方程左端的未来时刻变量值。但对于“隐式”及“半隐式”格式只是简单提及其概念和特点, 比如隐式格式为用未来时刻变量值求出未来时刻变量值, 它具有计算稳定、但计算复杂的特点。

在《数值天气预报》课程中, 需讲解大量的公式推导和讲解, 若不能配合简单而又形象的举例和图形, 学生尤其是本科生作为授课对象, 将很难理解和接受本课程中的相关内容, 讲课的效果也将大打折扣。在非线性不稳定计算的举例中, 教材中[2]虽给出了采用不同的初值和不同的差分方案对计算稳定性的影响, 但并没有清楚地列出其求解过程, 因此学生很难了解隐式格式差分的具体求解过程, 对教材中列出的“显示”和“隐式”格式的各自优缺点更是难以理解。因此, 需要对两种格式的计算过程进行相应的讲解, 尤其是隐式格式。此外, 教材中[2]对同一个微分方程构造的两个不同的差分方程中, 除隐式格式和显式格式的差异外, 还存在着对采用了不同的差分格式, 即显式格式采用中央差格式, 隐式格式采用前差格式。本课程[2]已清楚的讲解到中央差格式虽具有较高的计算精度, 但在时间差分计算时存在计算解的问题, 若初值取得不当, 则计算解会有较大的振幅。因此, 从逻辑上讲, 教材中给出的不同的差分方案的影响, 实际上不仅仅来源于显式格式和隐式格式的差异, 还来源于对时间微分采用不同差分格式的差异。这又加大了学生对显式格式和隐式格式特点的理解难度。

针对上述问题, 本文将以简单的一维非线性平流方程为例, 给出隐式格式差分方程的具体求解过程, 重新探讨非线性计算不稳定现象, 目的是使学生更好地了解显式格式和隐式格式差分方程的求解过程, 深刻理解两种格式各自的优缺点。

二、非线性计算不稳定的计算实例

以上两式与教材[2]基本一致, 不同的是这里的x取值范围并不到1。在大气科学中, 方程或模式的计算可在全球或某一纬圈上进行。在该情况下, 没有纬向侧边界条件。对于上式而言, 可认为u在x=1的取值等于u在x=0的取值, 也即循环边界条件。在构造上述微分方程相应的差分方程过程中, 对 (1) 式和 (2) 式分别采用显式格式和隐式格式:

其中上标n为第n步, 下标i为第i个格点, 。可见, 与教材中不同的是, (1) 式和 (2) 式中均取了前差格式, 这样可避免由于三个时间层计算而出现的计算解问题, 有利于问题的讨论更加集中。

同样给定两种不同的初值, 两者仅相差一个常数:

三、隐式格式的求解

显式差分方程 (3) 的求解过程即是将已知的n时刻u值代入等式右端算出等式左端未知的n+1时刻u值, 可见, 求解过程简单。至于隐式差分方程 (4) , 其求解过程, 较复杂。首先将 (4) 式写在[0, 1) 的。可见, a、b、c均为已知的第n步值。采用循环边界条件可得三元二次方程组:

将 (7) 式中的三式相加可得:X+Y+Z=2a+2b+2c, 再令m (2a+2b+2c) =d, 该d值也是已知的第n步值, (7) 式可化为三元一次方程组:

最终可利用已知的a、b、c和d值分别求得n+1步未知的Z、Y、X值:

再分别将其减去c、b和a值, 可得n+1步的u2n+1、u1n+1和u0n+1。由此可见, 隐式差分方程的求解过程较为复杂。需指出的是, 本文在[0, 1) 仅选取了3个格点, 若选取教材中的10个点 (Δx=0.1) , 则需在10个格点上写出10个差分方程, 并进行联立, 求解十元一次方程组, 其求解过程更为复杂。

四、计算结果及分析

图1a和1b分别给出了两个初值、两种计算方案的计算结果。初值取 (5) 式用显式方案 (3) 式的计算结果表明 (图1a实线) , 动能逐步增大, 在500步以后突然急剧增加, 出现按指数增加的趋势;但若给初值加上一个常数后 (图1b中实线) , 总动能在4m2/s2左右变化, 表明计算结果稳定。至于隐式方案 (4) 式, 无论取哪种初值, 结果均稳定。

至于产生如图1a中的不稳定现象, 仍可利用混淆误差理论进行解释, 即网格系不能正确分辨短波长的波动而导致不稳定。本文例子在[0, 1) 的一个周期范围内仅取三个格点, 采用循环边界条件, 即将[0, 1) 进行I=3等分。可见, 该网格系只能正确识别平均值0波、波长为的3/2波和波长2Δx为3Δx的1波波动。若波数k1=k2=3/2的两个波动相互作用, 则可产生0波和3波的波动。其中3波波动超出该网格系的识累, 从而可导致不稳定现象。该过程也可通过 (10) 式得到验证:

因此, 虽然本文的计算结果与教材中基本一致, 但举例十分简单, 这有利于学生的理解和接受。

五、结束语

本文通过简单的计算实例重新探讨了差分格式对非线性计算稳定性的影响。这里的“简单”, 主要指将一维非线性平流方程的时间偏导项统一地取成前差格式, 同时, 差分方程仅写在三个格点上。从而, 隐式格式差分方程的求解过程便成为三元一次方程组的求解过程。该求解过程比显式格式差分方程复杂, 但计算结果稳定, 充分体现出隐式格式和显式格式的优缺点。虽最终的计算结果与教材中[2]基本一致, 但本文举例更为简单、易懂, 且给出了详细的求解过程, 有助于学生自己动手推导和计算求解, 以加深其对显式和隐式格式的理解, 并深刻体会各自的优缺点。

摘要：数值天气预报是大气科学的一个重要分支, 是一门实用性很强的应用基础学科。数值天气预报本质是用数值方法求解非线性的大气运动方程组。但当采用格点差分来表示微分方程中的非线性项时, 易产生非线性计算不稳定现象。本文以一个简单的一维非线性平流方程的数值求解过程为例, 给出隐式格式差分方程的基本计算方法并重新演示了非线性计算不稳定现象。由于举例更为简单, 可加深学生对显式格式和隐式格式优缺点的理解。

关键词：大气科学,非线性计算不稳定,隐式格式

参考文献

[1]纪立人.数值天气预报发展进程中若干亮点的回顾及其启迪[J].气象科技进展.2011, 1 (1) :40-42.

[2]沈桐立, 田永祥, 葛孝贞, 陆维松, 陈德辉.数值天气预报[M].北京:气象出版社, 2003.9.

非线性计算不稳定篇2

问题描述及预备知识:

考虑不确定非线性切换系统

其中σ:[0, +∞) →I={1, 2LN}为切换规则, 是右连续的分段常值函数。x∈Rn是状态, d∈Rm为扰动输入, f i (⋅) , pi (⋅) 是已知的光滑函数, gi (⋅) 光滑且非奇异, δi (⋅) 是未知的光滑函数且满足:

ei是正实数, i∈I。

引理1考虑非线性切换系统

其中:x∈Rn是状态, d∈Rm为扰动输入, fi:Rn×Rm→R是关于x的局部Lipchitz函数, σ:[0, +∞) →I={1, 2LN}是右连续的分段常值函数.如果存在一组连续可微的正定函数iV (x) , i∈I和K∞类函数α1, α2与γ, 并且存在常数µ≥1, λ>0, 对∀x∈Rn, d∈Rm和∀p, q∈I满足:

若切换规则满足平均驻留时间, 则系统 (3) 关于d是输入对状态稳定的.

主要结果:

定理1对于系统 (1) , 存在一组连续可微的正定函数iV (x) , i∈I, 和K∞类函数α1, α2与ρ, 以及常数µ≥1, λ>0, 对∀x∈Rn, d∈Rm和∀p, q∈I满足:

则存在状态反馈控制器

在满足平均驻留时间的切换规则下, 闭环系统 (1) 关于d是输入对状态稳定的。

将iu的表达式代入上式得:

当时, 对于系统 (1) 的每一个子系统在

取并且切换规则满足平均驻留时间, 因此由引理1可得此结论成立。

摘要：本文针对一类含有结构不确定的非线性切换系统基于平均驻留时间的方法设计状态反馈控制器使得闭环系统是输入对状态稳定的。

关键词：非线性切换系统,不确定性,平均驻留时间,输入对状态稳定

参考文献

[1]张霞, 高岩, 夏尊铨.切换线性系统稳定性研究进展.控制与决策.Vol.25, No.10, Oct.2010.

[2]L.Vu*, D.Chatterjee, D.Liberzon.Input-to-state stability ofswitched systems and switching adaptive control.Automatica 43 (2007) 639-646.

[3]向峥嵘, 向伟铭, 陈庆伟.一类含扰动的非线性切换系统稳定性分析.控制与决策.Vol.23, No.1, 2008.

非线性计算不稳定篇3

关键词：干旱区,蒸发,有限元,排水系统

在蒸发强烈的地区分析排水问题时, 必须要考虑蒸发对地下水运动的影响[1]。但在排水计算模型中, 往往有许多被研究区域的物理量, 常常呈现为一复杂的函数, 当边界也较复杂时, 就很难得到该方程的解析解。而有限单元法则通过把研究域离散成有限个小单元, 用单元函数逼近总体函数, 从而把微分方程求解的问题化为一个求极值的问题, 最终得到研究区域研究对象的结果。对不稳定流有限单元法排水计算, 许多学者做了研究, 如王君连[2]推导出了补给项 (蒸发或灌溉) 为常数、有已知流量边界的二维不稳定流有限元方程, 李定方[3]推导出了多重补给条件的二维地下水稳定流有限元方程。本文应用Galerkin法推导了地下水蒸发与埋深呈线性关系时的二维不稳定流排水的有限单元方程, 并通过实例验证了方程的正确性。在此基础上, 用该有限元方程, 分析了新疆灌区排水系统的排水效果。

1 模型建立及求解

考虑蒸发补给时, 地下水的二维不稳定运动基本方程式可写为:

$\frac{\partial Η}{\partial t} = \frac{1}{μ} [\frac{\partial}{\partial x} (k_{x} Η \frac{\partial Η}{\partial x}) + \frac{\partial}{\partial y} (k_{y} Η \frac{\partial Η}{\partial y})] - \frac{ε}{μ} (1)$

式中:Kx, Ky为含水层x方向和y方向的平均导水率, m/d;μ为含水层平均给水度;H为以xoy平面为基准面的地下水位, m;t为排水历时, d;ε为地下水蒸发强度, m/d。

地下水蒸发强度ε与地下水埋藏深度Δ的关系一般用下式表示:

$ε = ε_{0} (1 - \frac{Δ}{Δ_{0}})^{n} (2)$

式中:ε0为地下水接近地表面时的蒸发强度, m/d, 其值接近水面蒸发强度;Δ为地下水埋深, m;Δ0为地下水蒸发为零时的埋深, m;n为蒸发指数, 一般为1～4。

取n=1, 因为Δ=T-H, T为自地面至不透水层的距离, m。

故:

$ε = ε_{0} - \frac{ε_{0}}{Δ_{0}} Τ + \frac{ε_{0}}{Δ_{0}} Η (3)$

将式 (3) 代入式 (1) 得:

$\frac{\partial Η}{\partial t} = \frac{1}{μ} [\frac{\partial}{\partial x} (Κ_{x} Η \frac{\partial Η}{\partial x}) + \frac{\partial}{\partial y} (Κ_{y} Η \frac{\partial Η}{\partial y})] - A Η - B (4)$

式中: $A = \frac{ε_{0}}{Δ_{0} μ} ‚ B = \frac{ε_{0}}{μ} - A Τ$ 。

将求解域剖分成四节点的矩形单元, 单元边长分别为2a和2b, 由Galerkin有限单元法的思想, 并认为渗流域总剩余量是单元剩余量之和[2], 则有下式成立:

$\begin{array}{l} \sum_{e = 1}^{Μ} \oint_{Γ_{e}} \frac{Η Ν_{L}}{μ} [k_{x} \frac{\partial \hat{Η}}{\partial x} \cos (n, x) + k_{y} \frac{\partial \hat{Η}}{\partial x} \cos (n, y)] d s - \\ \sum_{e = 1}^{Μ} \int_{e} [\frac{k_{x} Η}{μ} \frac{\partial \hat{Η}}{\partial x} \frac{\partial Ν_{L}}{\partial x} + \frac{k_{y} Η}{μ} \frac{\partial \hat{Η}}{\partial y} \frac{\partial Ν_{L}}{\partial y}] d x d y - \\ \sum_{e = 1}^{Μ} \int_{e} A \hat{Η} Ν_{L} (x, y) d x d y - \sum_{e = 1}^{Μ} \int_{e} B Ν_{L} (x, y) d x d y - \\ \sum_{e = 1}^{Μ} \int_{e} \frac{\partial \hat{Η}}{\partial t} Ν_{L} (x, y) d x d y = 0 (5) \end{array}$

式中: $\hat{Η} (x, y) = Η_{i} Ν_{i} (x, y) + Η_{j} Ν_{j} (x, y) + Η_{m} Ν_{m} (x, y) + Η_{n} Ν_{n} (x, y) = \sum_{L = i}^{n} Η_{L} (t) Ν_{L} (x, y)$ 为由Galerkin法构造的近似解表达式; $\oint_{Γ_{e}} \frac{Η Ν_{L}}{μ} [k_{x} \frac{\partial \hat{Η}}{\partial x} \cos (n, x) + k_{y} \frac{\partial \hat{Η}}{\partial x} \cos (n, y)] d s$ 是应用格林 (Greens) 公式, 把对面积的积分化为沿该面积为Γe周边的积分, 其含义是指沿有流量q流入的单元边界的积分;M为有限单元数, NL为形函数, HL为节点水头, m。

略去有限单元法的推导过程, 最后得到如下总方程:

$[Τ]^{(e)} [Η]^{(e)} = [F]^{(e)} (6)$

式中:

$\begin{array}{l} [Τ]^{(e)} = - \frac{Η}{3 μ} [\begin{matrix} (\frac{b}{a} k_{x} + \frac{a}{b} k_{y}) & (\frac{- b}{a} k_{x} + \frac{a}{2 b} k_{y}) & (\frac{- b}{2 a} k_{x} - \frac{a}{2 b} k_{y}) & (\frac{b}{2 a} k_{x} - \frac{a}{b} k_{y}) \\ (\frac{- b}{a} k_{x} + \frac{a}{2 b} k_{y}) & (\frac{b}{a} k_{x} + \frac{a}{b} k_{y}) & (\frac{b}{2 a} k_{x} - \frac{a}{b} k_{y}) & (\frac{- b}{2 a} k_{x} - \frac{a}{2 b} k_{y}) \\ (\frac{- b}{2 a} k_{x} - \frac{a}{2 b} k_{y}) & (\frac{b}{2 a} k_{x} - \frac{a}{b} k_{y}) & (\frac{b}{a} k_{x} + \frac{a}{b} k_{y}) & (\frac{- b}{a} k_{x} + \frac{a}{2 b} k_{y}) \\ (\frac{b}{2 a} k_{x} - \frac{a}{b} k_{y}) & (\frac{- b}{2 a} k_{x} - \frac{a}{2 b} k_{y}) & (\frac{- b}{a} k_{x} + \frac{a}{2 b} k_{y}) & (\frac{b}{a} k_{x} + \frac{a}{b} k_{y}) \end{matrix}] - \\ \frac{A}{9} [\begin{matrix} 4 a b & 2 a b & a b & 2 a b \\ 2 a b & 4 a b & 2 a b & a b \\ a b & 2 a b & 4 a b & 2 a b \\ 2 a b & a b & 2 a b & 4 a b \end{matrix}] - \frac{1}{9 Δ t} [\begin{matrix} 4 a b & 2 a b & a b & 2 a b \\ 2 a b & 4 a b & 2 a b & a b \\ a b & 2 a b & 4 a b & 2 a b \\ 2 a b & a b & 2 a b & 4 a b \end{matrix}] (7) \end{array}$

$\begin{array}{l} [Η]^{(e)} = [\begin{array}{l} Η_{i} \\ Η_{j} \\ Η_{m} \\ Η_{n} \end{array}] (8) \\ [F]^{(e)} = \frac{q l}{2 μ} [\begin{array}{l} 0 \\ 0 \\ 1 \\ 1 \end{array}] - \frac{B S}{4} [\begin{array}{l} 1 \\ 1 \\ 1 \\ 1 \end{array}] + \frac{1}{9 Δ t} [\begin{matrix} 4 a b & 2 a b & a b & 2 a b \\ 2 a b & 4 a b & 2 a b & a b \\ a b & 2 a b & 4 a b & 2 a b \\ 2 a b & a b & 2 a b & 4 a b \end{matrix}] [\begin{array}{l} Η_{i}^{-} \\ Η_{j}^{-} \\ Η_{m}^{-} \\ Η_{n}^{-} \end{array}] (9) \end{array}$

式中:S为矩形单元面积;H-i, H-j, H-m, H-n表示前期节点水头, 为已知值, 其余符号意义同前。

如区域的边界条件和初始条件给出, 则可按规则形成总体矩阵后, 利用Seidel迭代法求解方程的解。由于在计算前, 单元水头[H] (e) 并不知道, 为此需要采用“双重迭代法”, 即事前假定一组单元节点水头H-i, H-j, H-m, H-n, 则 $Η^{(e)} = \frac{1}{4} (Η_{i}^{-} + Η_{j}^{-} + Η_{m}^{-} + Η_{n}^{-}) ‚ Η^{(e)}$ 成为已知值。如果第一次计算结果不满足允许误差要求, 则H (e) 改用刚算出的一组节点水头取平均, 重新计算, 直到满足误差要求为止。

2 模型验证

为验证模型的正确性, 用该模型求解了文献[3]的排水问题。

已知某排水系统, 地面至不透水底层之距T=10 m, 地表面潜水蒸发强度ε0=0.006 m/d, 潜水蒸发为零值的地下水埋深Δ0=2.0 m, 土壤给水度μ=0.07, 渗透系数k=1.74 m/d, 水文地质特征值a2=183 m2/d, 排水沟长度L=S=120 m。计算1～5 d农田中点地下水位。

初始条件:H (x, y, t) =10 0<x<120, 0<y<120

边界条件:H (x, 0, t) =H (x, 120, t) =8

H (0, y, t) =H (120, y, t) =8, t>0

在划分计算区域时, 由于边界对称, 因此仅研究1/4个区域。取矩形单元边长2a=2b=5 m, 共剖分169个节点, 144个单元。时间步长Δt取等步长, 分别为Δt1=0.1 d和Δt2=0.05 d, 允许误差取0.05, 用VB语言将计算过程编制成程序进行计算。

根据已知条件, 把文献[3]的解析解和用本文中式 (6) 计算的关于农田中点处1～5 d的地下水位变化过程绘成图1。

由图1可知两种计算结果吻合较好, 因此, 可用本文的数学模型及程序计算蒸发影响的平面二维不稳定流排水计算问题。另外, 从图1可以看出, 时间步长分别为Δt1=0.1 d和Δt2=0.05 d时, 计算结果非常接近, 所以, 在计算时可取时间步长为0.1 d。

3 模型应用

西北干旱内陆区, 由于干燥少雨, 蒸发强烈, 当地下水位长期过高时, 就会形成土壤盐碱化, 因此, 该地区的排水系统主要是以降低地下水位为主。同时, 与南方圩垸灌区不同的是, 该地区的灌溉系统和排水系统是两个作用完全不同的系统。在有排水要求的灌区, 总是将灌溉农渠建在条田高处, 而将排水农沟建在条田低处, 两条农沟的间距一般就是条田的宽度。由于降雨稀少, 在淋洗土壤盐分时, 希望更多的水量先渗入土壤, 然后再进入排水系统, 以达到更好得冲盐和洗盐的目的。这样, 在外界蒸发影响的条件下, 如何确定合理的排水沟间距, 使之能在要求的时间内将地下水位降低到地下水临界深度将是灌区排水系统规划的重要问题。

本文以新疆巴州和硕县试验农田的排水系统为例, 应用推导的二维不稳定流有限单元方程计算了农田地下水的变化过程。

已知试验农田的排水沟长度L=240 m, 宽度S=100 m, 排水沟深2.5 m。试验农田纵向平均地面坡度为1/150 (图2中沿x轴方向) , 横向平均地面坡度为1/2 000, 由于横向坡度很小, 计算时取横向水平。为观测地下水位变化过程, 沿横向在试验农田中部每隔25 m设一个观测井。试验农田水系及观测井布置见图2。试验农田灌水时间为7月5日, 灌水后, 于每天上午10:00观测一次地下水位, 连续观测30 d。

排水计算时所需基本参数的确定如下所述:

(1) 根据钻孔资料表明, 地面至不透水底层之距T=15 m;

(2) 根据20 cm蒸发皿的观测资料, 测得7月5日～7月15日水面平均蒸发强度ε0=0.01 m/d;

(3) 土壤给水度μ和潜水蒸发为零值的地下水埋深Δ0根据文献[3]所述方法, 选取农田中部地下水面基本保持水平变化的数据进行计算。经过计算, Δ0=2.2 m, 土壤给水度μ=0.12;

(4) 渗透系数k按下式计算:

$k = 2 d_{10} e^{2} (10)$

式中:K为饱和导水率, cm/min;d10为有效粒径, mm;e为土壤孔隙比。

经过计算, K=0.225 cm/min。

计算初始条件:H (x, y, t) =16.6-x/150, 0<x<240, 0<y<100;

计算边界条件:H (x, 0, t) =H (x, 100, t) =14.1-x/150

$\frac{\partial Η (0, y, t)}{\partial x} = 0 ‚ Η (240, y, t) = 12.5, t > 0$

在划分计算区域时, 由于边界沿平行x轴对称, 因此仅研究1/2个区域。取矩形单元边长2a=15 m, 2b=5 m, 则共剖分187个节点, 160个单元。计算结果见图3～图5。

通过对图3～图5分析认为:

(1) 农田中部的地下水位 (y=50 m处) 是该农田排水系统的控制水位 (见图4) 。从图中可以看出, 地下水位在排水系统和地面坡度作用下, 其降落速度经历了由快到慢的过程, 前两天的累积降深平均为1.2 m, 后3天的累积降深为平均为0.9 m。到了第5天, 地下水位已基本降到了蒸发为零值的地下水位附近, 表明蒸发已经对地下水位的影响非常小了。另外, 从图中还可看出, 下游侧的排水沟对地下水位有显著影响的距离约为60 m。

(2) 地下水位的变化过程见图5。从图中可以看出, 在地面倾斜条件下, 地下水位在重力和排水沟共同作用下逐渐降低, 而且距离排水沟越近, 地下水位变化越大;排水1 d后, 农田最高地下水位为15.91 m, 潜水距离地面的最近距离约0.45 m, 排水5 d后, 农田最高地下水位降为14.16 m, 潜水距离地面的最近距离达到约2.14 m。通过分析认为, 该排水系统能够在5 d内将地下水位降低到蒸发影响很小的程度, 表明试验农田排水系统的排水效果较好。

(3) 试验农田的排水农沟间距为100 m时, 排水5 d后, 地下水位就可降低到2 m以下, 这对淋洗土壤盐分和降低因蒸发作用而使土壤表层积盐可起到积极作用。但是, 由于排水沟占地较多, 在新疆部分灌区, 对排水系统规划和施工时, 或者将排水沟间距放大, 或者减少排水沟级数, 只开挖干沟或支沟, 而使排水效果减弱。日积月累, 土壤便会次生盐碱化, 因此, 为有效治理土壤盐碱化, 应完善排水系统并选取合适的排水农沟间距。

参考文献

[1]王文焰, 李智录, 沈冰.对考虑蒸发影响下农田排水沟 (管) 间距计算的探讨[J].水利学报, 1992, (7) :23-28.

[2]王君连.工程地下水计算[M].北京:中国水利水电出版社, 2004.

[3]沙金煊.农田不稳定排水理论与计算[M].北京:水利水电出版社, 2004.

[4]马英杰.地面倾斜条件下考虑蒸发影响的农田不稳定流排水计算[J].中国农村水利水电, 2008, (1) :56-62.

[5]张蔚臻.地下水非稳定流计算及地下水资源评价[M].北京:科学出版社, 1983.

[6]何光渝.VB常用算法大全[M].西安:西安电子科济大学出版社, 2001.

[7]张蔚臻论文集编写组.张蔚臻论文集[M].武汉:武汉大学出版社, 2002.

[8]Benoit G R.The effect of soil surface condition on evaporation ofsoil water[J].Soil Sci.Amer.Proc, 7:495-498, 1963.

[9]Brutsaert W.Research note onthe first and secondlinearation of theBoussinesq equation[J].Geophysics, 11:549-554, 1965.

[10]Cooley R L.Some new procedures for numerical solution of vari-ably saturated flow problems[J].W, R, R19 (5) , 1983.

【非线性计算不稳定】推荐阅读：