CN105224806B - 一种获取非晶态物质径向分布的方法 - Google Patents

Abstract

本发明公开了一种获取非晶态物质径向分布的方法，属于材料科学领域。其步骤如下(1)确定所需要获得的径向分布的函数，给定初始近似x₀和允许误差ε；(2)计算非线性方程组f(x_k)和▽f(x_k)，如果有||f(x_k)||<ε或者||▽f(x_k)||<ε，则停止计算，x_k为满足精度的解，否则转向(3)；(3)确定非晶态物质径向分布解的迭代序列{x_k}；(4)计算下降方向p_k，Hessian矩阵函数s_k，q_k，如果转向(2)，否则取H_k+1＝I，单位矩阵I，即H_k+1＝I，然后转向(2)；它可以避免计算时候的运算中断，正确性高、效率高。

Description

一种获取非晶态物质径向分布的方法

技术领域

本发明涉及材料科学领域，更具体地说，涉及一种获取非晶态物质径向分布的方法。

背景技术

我们在进行材料科学研究时，往往需要对物质的径向分布进行研究，非晶态物质的径向分布比较复杂，如硅等，N元非晶态物质的径向分布函数为：

其中偏径向分布函数

RDF_αβ(r)＝4πr²ρ_αβ(r) (2)

式中ρ_αβ(r)指以第α种原子为原点在距离r处的第β种原子的数密度。显然要求RDF_α只需要求得沿径向的原子数密度函数ρ_αβ(r)即可。之前人们是通过傅里叶变换来得到原子数密度的，但是这种方法必须对未知衍射强度做出假设，一方面不够严谨，另一方面得到结果有失真。为此人们又利用最大熵原理通过求解非线性方程组得到径向分布函数。该非线性方程组为：

其中

r_j＝jΔr为第j个球壳的半径，Δr为球壳厚度，ρ_β,0为非晶态物质样品内第β种原子的平均原子数密度，δ_αβ为Kronecker符号，n是参与散射的原子个数，k_i为第i个散射矢量，I(k_i)是散射强度，f原子散射因子，L为实验数据点个数。一般地J≤200；ρ_β,0，k_i，I(k_i)，f均为实验所测，n根据散射曲线拟合求得。具体拟合方法可参考B.E.Warren,X-RayDiffraction,chapter 10，Dover Publication,New York,1990。

我们的目标就是求解方程组(3a)和(3b)，从而得出μ_αβ和μ_i，将其代入(4)式，进而求得ρ_αβ(r)。然而，通常情况下数据点个数L为10³的量级，因此未知数多、数据量庞大，再加上方程为指数形式的非线性方程，求解十分困难。

变尺度法是求解非线性方程组的重要而极为有效的方法之一，由于该方法避免了计算二阶导数矩阵及对其求逆的过程，因此计算量大幅减小；其比梯度法和共轭梯度法的收敛范围更广，收敛速度更快，在处理高维最优化问题有明显优势。在诸多变尺度方法中，目前公认的最成功的为BFGS子算法。下面我们简略介绍之，指出BFGS算法在实际计算时所存在的问题。

对多维问题我们用向量来表示：

其中向量函数f(x)的任意一个分量f_i(x)都为向量x的非线性的函数，那么非线性方程组可写为：

f(x)＝0 (6)

我们称f(x)为定义在上且取值于Rⁿ的向量值函数，将它简记为f：若存在满足方程组(2)，则称x＊为方程组的解。

在无法求出解析解的情况下，只能数值解之。然而，必须明白的是方程组的解并不总是能准确找到，这一方面取决于方程组的解是否存在，另一方面取决于采用的算法和收敛精度，以及初始值的选取。如果在计算过程中有序列{x_k}使下式成立：

即在Rⁿ的2-范数意义下有：

||x_k-x*||→0(k→∞) (8)

则称迭代序列{x_k}收敛于x＊，即x_k(k→∞)为方程组的解。实际求解中，这是自相矛盾的，因为如果知道了x＊，那就没有必要再计算(7)或者(8)式了。那么如何判断序列{x_k}是否正在向x＊收敛呢？实际上是通过f(x)的范数来判断的，如果它收敛于某一恒量，即：

||f(x_k)||<ε (9)

ε为允许误差，那么此时的x_k便为方程组(6)的解。

如上所述，迭代计算的目的是让目标函数快速、稳定地收敛于ε。变尺度法中的迭代序列{x_k}由下式给出：

w_k为步长因子，p_k为下降方向，H为Hessian矩阵。在BFGS算法中，H按下式构造：

其中：

而初始Hessian矩阵H₀为单位矩阵I。Hessian矩阵的对称正定性要求下式必须被满足：

尽管在理论上(13)式总是成立的，但是在实际计算过程中计算机的存储精度是有限的，因此当为一极小量时，在存储时就变成了甚至在中间计算过程类似的问题导致此时(13)式就不再成立，Hessian矩阵将不再是正定矩阵，运算要么被中断要么给出错误的结果，进而无法有效求解方程组(6)。现有公开中也没有合适的方法解决这个问题。

发明内容

1.要解决的技术问题

为了克服现有技术中获得非晶态物质径向分布困难、计算易出错、计算量大、效率低的问题，本发明提供了一种获取非晶态物质径向分布的方法。它可以避免计算时候的运算中断，正确性高、效率高。

2.技术方案

本发明的目的通过以下技术方案实现。

一种获取非晶态物质径向分布的方法，步骤如下：

(1)确定所需要获得的非晶态物质径向分布的函数，给定初始近似x₀和允许误差ε；

(2)计算f(x_k)和▽f(x_k)，如果有||f(x_k)||<ε或者||▽f(x_k)||<ε，则停止计算，x_k为满足精度的解，否则转向(3)；

(3)确定非晶态物质径向分布解的迭代序列{x_k}；

w_k为步长因子，p_k为下降方向，H为Hessian矩阵；

H按下式构造：

其中：

(4)如果直接转向(2)，否则取H_k+1＝I，I为单位矩阵，然后转向(2)。

由此获得非晶态物质径向分布函数的解x_k。

更进一步，步骤(1)为通过最大熵原理求解非线性方程组得到径向分布函数的非线性方程组为：

其中

μ_αβ,μ₁,μ₂,…为所求的未知数，r_j＝jΔr为第j个球壳的半径，Δr为球壳厚度，ρ_β,0为非晶态物质样品内第β种原子的平均原子数密度，δ_αβ为Kronecker符号，n是参与散射的原子个数，k_i为第i个散射矢量，I(k_i)是散射强度，f原子散射因子，L为实验数据点个数，J≤200。

更进一步的，该方法可用于求解任何大规模的非线性方程组。

3.有益效果

相比于现有技术，本发明的优点在于：

(1)本发明将BFGS算法和梯度算法有机结合起来，强有力地保障了Hessian矩阵的正定性，增强了算法的数值稳定性，在整个求解运算过程中不会因为计算机的存储精度而被迫中断或者给出错误的结果；

(2)使用变尺度算法针对非晶体径向分布的获得，提高了运算效率，大幅缩短计算耗时、降低计算成本。

附图说明

图1为本发明的变尺度算法流程图；

图2为利用本发明方法得到的非晶硅径向原子数密度。

具体实施方式

下面结合说明书附图，对本发明作详细描述。

现上述技术背景中所提及的不能被满足的问题，本发明在算法上进行了改进，避免了在(13)式无法被满足的情况下运算的中断或者给出错误结果。在运算过程中使用Hessian矩阵求解，当Hessian矩阵为非正定矩阵时，将其重置为单位矩阵。

①迭代计算过程中每一步的x_k的中间值被存储，检验由此计算出来的是否满足(12)式，如满足则进行下一步迭代，如不满足则转入①；

②将Hessian矩阵的中间矩阵H_k重置为单位矩阵I，即H_k＝I，由此计算出x_k+1进行下一步迭代。

在②中进行了Hessian矩阵重置为单位矩阵后，紧接着进行的一步迭代实际上为梯度法，即由BFGS算法转换为梯度法算法，这一转换避免了由于计算机存储精度导致的运算中断或者错误，并能在下一步迭代再次转入BFGS算法。由此可在一定程度上确定非晶态物质的微观结构。

如图1所示：

(1)N元非晶态物质的径向分布函数为：

其中偏径向分布函数

RDF_αβ(r)＝4πr²ρ_αβ(r) (2)

式中ρ_αβ(r)指以第α种原子为原点在距离r处的第β种原子的数密度非线性方程组为：

其中

r_j＝jΔr为第j个球壳的半径，Δr为球壳厚度，ρ_β,0为非晶态物质样品内第β种原子的平均原子数密度，δ_αβ为Kronecker符号，n是参与散射的原子个数，k_i为第i个散射矢量，I(k_i)是散射强度，f原子散射因子，L为实验数据点个数。一般地J≤200；ρ_β,0，k_i，I(k_i)，f均为实验所测，n根据散射曲线拟合求得。

给定初始近似x₀和允许误差ε。

求解时向量表示如下：

其中向量函数f(x)的任意一个分量f_i(x)都为向量x的非线性的函数，那么非线性方程组可写为：

f(x)＝0 (6)

我们称f(x)为定义在上且取值于Rⁿ的向量值函数，将它简记为f：若存在满足方程组(2)，则称x＊为方程组的解。

在无法求出解析解的情况下，只能数值解之。然而，必须明白的是方程组的解并不总是能准确找到，这一方面取决于方程组的解是否存在，另一方面取决于采用的算法和收敛精度，以及初始值的选取。如果在计算过程中有序列{x_k}使下式成立：

即在Rⁿ的2-范数意义下有：

||x_k-x*||→0(k→∞) (8)

则称迭代序列{x_k}收敛于x＊，即x_k(k→∞)为方程组的解。实际求解中，这是自相矛盾的，因为如果知道了x＊，那就没有必要再计算(7)或者(8)式了。那么如何判断序列{x_k}是否正在向x＊收敛呢？实际上是通过f(x)的范数来判断的，如果它收敛于某一恒量，即：

||f(x_k)||<ε (9)

ε为允许误差，那么此时的x_k便为方程组(6)的解。

(2)计算f(x_k)和▽f(x_k)，如果有||f(x_k)||<ε或者||▽f(x_k)||<ε，则停止计算，x_k为满足精度的解，否则转向(3)。

(3)确定非晶态物质径向分布解，计算时候的迭代序列{x_k}由下式给出：

w_k为步长因子，p_k为下降方向，H为Hessian矩阵。在BFGS算法中，H按下式构造：

其中：

计算出第k次近似x_k。

(4)计算p_k，s_k，q_k。如果直接转向(2)计算，否则取H_k+1＝I，单位矩阵I，即H_k+1＝I，然后转向(2)计算。

由此获得结果非晶态物质径向分布函数的解x_k。

具体的实施方式如下：

对于非晶硅，我们在实验上测得了散射强度曲线，平均原子数密度初始近似取x₀＝(2.0119,0,0,…,0),允许误差ε＝0.1.按照上述方法成功求解得出了如图2所示的径向原子数密度。

以上示意性地对本发明创造及其实施方式进行了描述，该描述没有限制性，在不背离本发明的精神或者基本特征的情况下，能够以其他的具体形式实现本发明。附图中所示的也只是本发明创造的实施方式之一，实际的结构并不局限于此，权利要求中的任何附图标记不应限制所涉及的权利要求。所以，如果本领域的普通技术人员受其启示，在不脱离本创造宗旨的情况下，不经创造性的设计出与该技术方案相似的结构方式及实施例，均应属于本专利的保护范围。此外，“包括”一词不排除其他元件或步骤，在元件前的“一个”一词不排除包括“多个”该元件。产品权利要求中陈述的多个元件也可以由一个元件通过软件或者硬件来实现。第一，第二等词语用来表示名称，而并不表示任何特定的顺序。

Claims (2)

1.一种获取非晶态物质径向分布的方法，步骤如下：

(1)确定所需要获得的非晶态物质径向分布的函数，给定初始近似x₀和允许误差ε；步骤(1)为通过最大熵原理求解非线性方程组得到径向分布函数的非线性方程组为：

其中

μ_αβ,μ₁,μ₂,…为所求的未知数，r_j＝jΔr为第j个球壳的半径，Δr为球壳厚度，ρ_β,0为非晶态物质样品内第β种原子的平均原子数密度，δ_αβ为Kronecker符号，n是参与散射的原子个数，k_i为第i个散射矢量，I(k_i)是散射强度，f原子散射因子，L为实验数据点个数，J为所要计算的球壳的个数，J≤200；

(2)计算向量函数f(x_k)和▽f(x_k)，如果有||f(x_k)||<ε或者||▽f(x_k)||<ε，则停止计算，x_k为满足精度的解，否则转向(3)；

(3)确定非晶态物质径向分布解的迭代序列{x_k}；

w_k为步长因子，p_k为下降方向，H为Hessian矩阵；

H按下式构造：

其中：

H₀＝I，I为单位矩阵；

(4)如果直接转向(2)，否则取H_k+1＝I，I为单位矩阵，然后转向(2)；

由此获得非晶态物质径向分布函数的解x_k。

2.根据权利要求1所述的一种获取非晶态物质径向分布的方法，其特征在于：该方法可用于求解任何大规模的非线性方程组。