细节：稳健的数值计算代码

1月前浏览192

作为一名数值算法开发人员，今天想聊一聊开发中的一些细节问题，我以C++代码为例，说一下在数值算法开发中需要注意的一些代码细节问题，这些细节问题直接决定你的代码稳健程度。

首先大家都知道，浮点数是有误差的，float的误差大概在6~9位有效数字，而double的误差大概在15~17位有效数字，由于这些误差的存在，就会造成一些数值上的问题，而有些不符合数学计算规则的异常数据的出现也会影响很多计算或判断流程。今天看看以下的情况：

任何数字与NaN比较均返回False

这段代码我们本来是想判断b是否在[-1，1]的区间内，然而由于sqrt（-1）返回NaN，所以最后我们会得出Hey this b is in this range的结果，这显然是不对，大概率这会导致其走上一条错误的路继续计算。一种更为成熟稳定的做法就是改变判断的方式，看如下的代码

那么同样是判断b是否在区间内，上述代码则会返回

Hey the b is not in this range
This b is: -nan

显然这是一种更为合理的方式。在有限元中编程、计算几何编程中，我们常常会有开根号、做除法的操作，当遇到异常数据的时候，第二种方法“将想要判断场景放在正向的判断中”往往可以更好的Debug，避免异常情况的进一步扩散。

非常危险的情况就是除0

假设我们要做如下计算

当与向量是垂直的时候，那么就是一种除零的情况，而这个计算在计算线段与平面是否相交的时候应用广泛，比较好的做法就是在计算除法之前先判断一下

double tmp = Dot(n, B-A);
if (IsSame(tmp, 0.0))
{
    if （onPlane(A, P))
    {
        reutrn true;
    }
    return false;
}
else
{
  if 0<=((d-Dot(n,A))/tmp)<=1
  {
      return true;
  }
}
// t not in [0, 1]
return false;

这里其实可以看到基本计算并没有改变，无非是增加了对于与是否垂直的判断，由于为平面P的法向，则只需要判断AB线段是否有一个点在这个平面上，就可以判断AB线段是否完全在这个平面上。另外当不垂直的情况下，则通过t是否在区间[0,1]内进行判断。

注意舍去误差

请看下边这个例子：

理论上b不为0，然而我们看一下输出的结果

大家千万不要忽视这个，觉得自己工作中不会遇到这样的问题，比如铁的弹性模量大概是 , 一根铁桁架假设那么 ,数量级已经很大了，在多次计算后可能就会出现上述的问题。

尽可能用相对误差

比如下面的代码if (Abs(x-y)<=epsilon那么一般情况下来说更好的都是以下的这种相对误差的写法if (Abs(x-y)<=epsilon*Max(x,y))

最后说一句，报错信息多点，总是好的！

来源：大狗子说数值模拟

结构仿真小知识：点对面及面对面

我工作后的第一任领导跟我说“魔鬼都在细节中”，这句话对我受益匪浅，所以这个系列我会说一说一连串的小知识点，这些大多数我们在仿真中会用到，但可能并不会关注，可有时候仿真的正确性就是由这些小的细节决定的今天说一下结构仿真中的接触，接触中有两种不同的接触计算方法，分别是点对面（nodetosurface）与面对面（surfacetosurface）。注意这里是两种不同的算法，或者准确的说是接触计算离散的两种不同算法，很多同志经常会把这个和设置接触的时候选择的对象混淆，觉得应该有点对线、线对线、线对线、点对面、面对线、面对面这么多种类型。事实上大多数情况，无论我们使用点对面算法还是面对面算法，我们通常都是选择两个几何面或者两个网格面，当然点对面顾名思义，他的算法天生支持其选一堆点作为从节点。下图为Abaqus中这两种方法的选择，在Interaction中选择完主从面后即可设置离散方法（Discretizationmethod）。那么现在我们就分别说一下点对面与面对面算法。首先点对面算法是最早出现也最成熟的算法，逻辑很简单，即从点向主面进行投影，然后得投影点到该从点的方向向量，通过这个方向向量与主面外法线，我们就可以得到从点与主面是处于分离状态，还是穿透状态，如果处于穿透状态，那么我们就要进行接触刚度的计算，针对每一个这样的点-面对进行依次的计算，用一行公式表达就是：,,其中为法向的间隙，当其小于0的时候代表穿透，真实世界不会发生，所以这是需要引入接触刚度，进行非线性求解，具体接触刚度怎么计算，这里先不详谈，后续会有专门讲计算接触力学的“大知识”专栏。从上述描述上其实我们可以看到一个特点，我们说从面的节点向主面投影，那么这里有两个问题：从面不能穿透主面，但好像没有任何一条公式表示主面不能穿透从面？从面如果很大，主面如果不够大，是否会存在从节点投影不到主面的情况？这两个问题的答案都是“是的”，主面理论上可以随意的穿透从面，而且像下边这个典型的案例，其实点对面接触将失效，因为虽然这两个面贴到一起了，但是由于都是单个单元，从面上的四个节点无法投影到主面上，也就无法人为接触是否发生了。这也是很多书上、帮助文档上都会有一条，“选择接触主面的时候建议选用大一些的面作为主面“。那么我们有没有办法规避这两个问题呢，对于问题1，诸如Ansys与Abaqus均提供了双向搜索的一种方式，Ansys中称作对称的接触（Symmetric）,而学术界则一般称其为双向接触（two-pass），理解起来也十分容易，那既然现在不保证主面不穿透从面，我就再把主从调换一下，再检测一下，再算一次，这回的原来的主面变成了从面也就不能穿透原来的从面了，这也就是双向的含义。而对于问题2，一个非常简单的方法就是，细化网格，把从面的网格细化了，这样从面的点多了，有些落在红色的主面投影范围内也就好了，这则是很多资料说“从面网格画细一点”的原因。另外一种做法是，诸如在Abaqus中他会针对这种情况去补充一些点，这里可以看到当选择了Nodetosurface算法后，他会出现Usesupplementarycontactpoints（使用补充的接触节点），一般来说我们默认的既可，那么Abaqus会自动再这种情况的地方，多做一些补充性的从节点，以避免这种漏检测的情况出现。而面对面算法，相对点对面就会复杂很多，相比点对面算法，这里则是将从面单元面投影到主面上，这个投影和面的选择又有很多玩法，留到未来计算接触力学中详谈，总之目前来说应用最为广泛的一种方法叫Mortar法（我并不知道这东西怎么翻译成中文），核心意思就是选择一个面作为基准面，主面和从面都向这个“中间”面进行投影，然后投影的主从面区域进行相互切割，并进行积分计算接触刚度、接触力等。面对面算法仍然符合点对面算法“无穿透”的基本数学公式，只是这里是通过一个积分来保证无穿透，简单的可以理解为从面单元上任意一点的穿透量（插值得来）在这主从两个面投影的交叉区域上的积分是“无穿透”的，说起来有点拗口，用个图来表示：点对面的感觉，红色点穿透了，穿透量是每个点到上面表面的投影面对面的感觉，穿透量是红色区域穿透量的积分，也就是红色和边界围出来那个小面积从上边这些中大家其实已经可以感觉到，面对面算法一定计算量更大一些，而且更准确一些，且不会出现点对面那种红、蓝主从检测不到的情况，因为使用面对面算法，那两个面是有接触面积的。也正是因为这些特性，面对面由于考虑了一个整体的效应，其接触面上的接触压力、应力分布都会光顺很多，但点对面（NTS）算法则会跳来跳去（即使你的面很光顺），对于二阶单元这个效应则会更加明显，这也是为何很多软件在点对面接触的时候并不推荐使用二阶单元。当然成熟的商业有限元分析软件或多或少都有考虑修正这些应力不光顺，比如Abaqus就提供了一些Smooth算法来消除这些效应，也同时提高收敛性，这个先提一下，以后可以单拎出来详谈，很有意思，还涉及到一些计算几何的知识。好的，今天就说到这，希望大家在以后的工作中，知道要用哪种接触，大多数情况用程序默认的（隐式：面对面，显式：点对面），少数情况为了计算效率、收敛性可以切换成点对面，但是也不要期盼在用点对面的时候可以多么正确的评估接触区域的接触压力与应力水平。来源：大狗子说数值模拟