在 COMSOL Multiphysics® 中曲线拟合解数据

2021年 7月 27日

在 COMSOL Multiphysics® 软件中求解一个模型后,我们可能希望将求解数据拟合到在仿真域中定义的一组函数中。在之前的博文中,我们解释了如何将离散试验数据拟合到曲线上。今天,我们将考虑对连续解数据进行拟合。之后,我们将介绍正交性的概念,并解释如何将解数据拟合到一组正交函数中,从而简化为简单方便的后处理操作。

曲线拟合的好处

在上一篇关于拟合离散试验数据的博客文章中,我们解释了如何使曲线拟合有助于随后的仿真工作。例如,如果我们尝试使用原始实验数据直接定义材料属性,那么数据中的统计波动会使求解器更难以收敛。此外,将离散数据拟合到平滑函数可以获得该函数的高阶导数,而尝试对原始数据进行数值微分可能会产生很大干扰且容易出错。

显示原始实验数据的图,可视化为黑点,拟合三次多项式,可视化为红色曲线
将原始实验数据(黑点)拟合为三次多项式(红曲线)。

对 COMSOL Multiphysics 模型中以前研究的解数据进行曲线拟合会怎么样?一个直接的好处是数据压缩。如果我们可以找到与完整解有良好一致性的函数的线性组合,那么我们可以通过共享几个函数系数的值来描述来自完整解的大部分信息。(请注意,对于数据压缩,没有放之四海而皆准的方法;在 COMSOL 知识库中,减少模型中存储的解数据量条目介绍了其他几种方法,例如使用探针表和选择)。

此外,对连续解数据进行曲线拟合可以方便地估计解的高阶空间导数。COMSOL® 软件中的许多物理场接口使用了由分段拉格朗日多项式组成的有限元离散化,默认情况下,这些多项式通常是二阶的。在这种情况下,导数不一定是连续的,三阶或高阶导数总是零。

最小二乘拟合概述

让我们首先看看简单最小二乘拟合的基本数学原理。在我们用 Ω 表示的某个模拟域(或边界或边)中,假设我们的 COMSOL Multiphysics 模型已经求解了场变量 u。我们的目标是通过将其视为一组预定义函数的线性组合来近似 u 在该域上的值,

u(\mathbf{r}) \approx \sum_{i=1}^N c_i f_i(\mathbf{r})

其中,fi 是一组已知函数,ci 是我们必须求解的系数。

确定这些未知系数的常用方法是,求 u 和它的近似值之差 L2 范数,

\int_\Omega \left(u(\mathbf{r})-\sum_{i=1}^N c_if_i(\mathbf{r})\right)^2\textrm{d}\Omega

然后,找到对应于这个 2 范数的局部最小值的i

在 Ω 内的局部极小值处,未知系数的导数必须等于零,

\frac{\partial}{\partial c_j}\left[\int_\Omega \left(u(\mathbf{r})-\sum_{i=1}^N c_if_i(\mathbf{r})\right)^2\textrm{d}\Omega\right]=0

其中,使用下标 j 以避免与求和的符号混淆。

假设积分号下的所有函数都可导,通过应用链式法则,它可以简化为

2\int_\Omega \left(u(\mathbf{r})-\sum_{i=1}^N c_if_i(\mathbf{r})\right)\left(-f_j(\mathbf{r})\right)\textrm{d}\Omega=0

约掉系数 2,交换积分和离散求和的顺序,得到

\sum_{i=1}^N c_i \int_\Omega f_i(\mathbf{r}) f_j(\mathbf{r}) \textrm{d}\Omega = \int_\Omega u(\mathbf{r})f_j(\mathbf{r})\textrm{d}\Omega

最后一个结果是一组有 N 个系数未知的 N 维线性方程组。例如,如果有三个多项式,那么我们可以把这个结果写成矩阵形式,

\begin{aligned}
&\begin{bmatrix}
A_{11} & A_{12} & A_{13}\\
A_{21} & A_{22} & A_{23}\\
A_{31} & A_{32} & A_{33}
\end{bmatrix}
\begin{bmatrix}
c_1\\
c_2\\
c_3
\end{bmatrix}
=
\begin{bmatrix}
B_1\\B_2\\B_3\\
\end{bmatrix}\\
&A_{ij} \equiv \int_\Omega f_i(\mathbf{r})f_j(\mathbf{r})\textrm{d}\Omega
\qquad B_{j} \equiv \int_\Omega u(\mathbf{r})f_j(\mathbf{r})\textrm{d}\Omega
\end{aligned}

 

因此,求解这个方程组的未知系数需要计算 N(N+3)/2 次积分(注意左边的矩阵是对称的,即 i f j 的积分等于 i f j 的积分)和一组N个线性代数方程的解。

正如我们将在下面的小节中看到的,如果我们定义相应的函数 i,使它们在我们试图拟合解数据的区域中是正交的,那么系数i的计算就会大大简化。

正交函数简介

在讨论正交性的概念时,我们必须注意不要给出一个过于狭隘的定义。一般而言,如果 \langle u,v\rangle = 0,属于一个向量空间的两个不同向量 uv 在内积下是正交的 ,内积的确切定义可以根据向量空间而改变。

由于我们讨论的是在 n 维实坐标空间(\mathbf{r}\in\mathbb{R}^n)的某个区域 Ω 上,一个空间变化的解向量对一组函数 f1(r), f2(r),…fN(r) 的曲线拟合,因此我们可以得到更具体的表达式,将这些函数的内积定义为

\langle f_i,f_j\rangle \equiv \int_\Omega f_i(\mathbf{r})f_j(\mathbf{r})w(\mathbf{r})\textrm{d}\mathbf{r}

其中,w(r) 为某个尚未定义的权函数,对于 Ω 中的所有 r,限制为:w(r) > 0。

例如,f_1(x) = 1,f_2(x) = \sin x ,和 f_3(x) = \cos x 在内积下正交。

\langle f_i, f_j \rangle \equiv \int_{-\pi}^{\pi}f_i(x)f_j(x)\textrm{d}x

也就是说,Ω 的区域是一维区间 x\in[-\pi,\pi] ,权函数就是w(x)=1。我们可以通过计算积分来证明这些函数的正交性,权函数就是 w(x)=1。我们可以通过计算积分来证明这些函数的正交性

\begin{aligned}
\int_{-\pi}^\pi (1)^2 \textrm{d}x &= 2\pi \qquad \int_{-\pi}^\pi (\sin x)^2 \textrm{d}x = \pi \qquad \int_{-\pi}^\pi (\cos x)^2 \textrm{d}x = \pi \\
\int_{-\pi}^\pi (1)(\sin x) \textrm{d}x &= 0 \qquad \int_{-\pi}^\pi (1)(\cos x) \textrm{d}x = 0 \qquad\int_{-\pi}^\pi (\sin x)(\cos x) \textrm{d}x = 0\\
\end{aligned}

事实上,我们可以进一步扩展这个概念,看到函数 sin(2x), cos(2x), sin(3x), cos(3x) 等在同一个内积下也是正交的。这些三角函数的正交性是傅立叶级数展开的基础。

正交函数的最小二乘拟合

现在,让我们利用我们对正交函数的了解,重新审视线性最小二乘拟合问题,正如我们之前所看到的,该问题简化为求解系数ci的集合,使得

\int_\Omega \sum_{i=1}^N c_i f_i(\mathbf{r}) f_j(\mathbf{r}) \textrm{d}\Omega = \int_\Omega u(\mathbf{r})f_j(\mathbf{r})\textrm{d}\Omega

现在,假设函数fi都是在以下内积下是正交的

\langle f_i, f_j \rangle \equiv \int_\Omega f_i(\mathbf{r})f_j(\mathbf{r})\textrm{d}\Omega

你可能已经注意到了一些技巧:在谈到 n 维向量空间的内积时,我们在内积的定义中使用了微分元 dΩ,但是 dΩ 将根据我们选择的坐标系(如笛卡尔坐标系、圆柱极坐标系、球极坐标系等)而采用不同的表达形式。正如我们将在后面的示例中看到的那样,为了最大程度的方便,我们可以选择一组在内积下正交的函数,其权函数等于所用坐标系的雅可比行列式。

由于函数的正交性,所有的非对角线项在 N×N 矩阵的左边消失了,现在 N 个方程组现在只是一组 N 个表达式的集合,可以直接得到 ci 的值

c_i = \frac{\int_\Omega u(\mathbf{r})f_i(\mathbf{r})\textrm{d}\Omega}{\int_\Omega f_i(\mathbf{r}
)^2\textrm{d}\Omega}

如果我们选择将这些函数规范化,使上述表达式的分母为 1,我们就会得到看起来更简单的表达式

c_i = \int_\Omega u(\mathbf{r})f_i(\mathbf{r})\textrm{d}\Omega

因此,将解数据拟合到一组 N 个函数的任务已经简化为计算 N 个不同的积分。在 COMSOL® 软件中,我们可以轻松定义积分 耦合,以计算在任何区域、边界、边或几何点集合的积分。(我们甚至不需要在定义积分 耦合后重新运行研究;只需单击“更新解” 即可。)

在接下来的内容里,我们将研究一个更实际的示例: 将镜子的变形拟合到一组 Zernike 多项式。

泽尼克多项式

光学中最常用的正交函数之一是泽尼克多项式,它是圆的径向坐标和平面角的函数,

Z_n^m(\rho,\theta) = N_n^m R_n^{|m|}(\rho)M(m\theta)

其中,

  • ρ 是径向坐标(0\leq\rho\leq 1)
  • θ 是方位角 (0\leq \theta < 2\pi)
  • N_n^m 是归一化项
  • R_n^{|m|}(\rho) 是径向项
  • M(m\theta) 是子午项或方位项
  • n 是径向指数 (n \in \left\\{0,1,2,\dots\right})
  • m 是子午线或方位角指数(对于给定的 n,m \in \left
    {-n, -n+2, -n+4, \dots , n-2, n\right}
     )

需要注意的是,对于如何定义泽尼克多项式以及如何解释其指数,有几种不同的标准格式。COMSOL Multiphysics 中使用的函数定义使用两个独立的指数来描述函数的径向和方位角相关性,符合 ANSI 和 ISO 标准(参考文献 1、2)。泽尼克多项式通常是针对单位圆定义的,但可以通过将 ρ 替换为 ρR来定义任何其他半径为 R 的圆(尽管这会影响归一化)。

下面列出了部分泽尼克多项式。

表达式 通用名
Z_0^0 1 平移
Z_1^{-1} 2\rho \sin(\theta) 垂直倾斜
Z_1^1 2\rho \cos(\theta) 水平倾斜
Z_2^{-2} \sqrt{6}\rho^2 \sin(2\theta) 斜散像差
Z_2^0 \sqrt{3}\left(2\rho^2-1\right) 散焦
Z_2^2 \sqrt{6}\rho^2 \cos(2\theta) 散光
Z_3^{-3} \sqrt{8}\rho^3 \sin(3\theta) 斜向三叶草像差
Z_3^{-1} \sqrt{8}\left(3\rho^3-2\rho\right)\sin(\theta) 垂直慧差
Z_3^1 \sqrt{8}\left(3\rho^3-2\rho\right)\cos(\theta) 水平慧差
Z_3^3 \sqrt{8}\rho^3 \cos(3\theta) 水平三叶草像差
Z_4^{-4} \sqrt{10}\rho^4\sin(4\theta) 斜四叶草像差
Z_4^{-2} \sqrt{10}\left(4\rho^4-3\rho^2\right)\sin(2\theta) 斜次阶散光像差
Z_4^0 \sqrt{5}\left(6\rho^4-6\rho^2+1\right) 球差
Z_4^2 \sqrt{10}\left(4\rho^4-3\rho^2\right)\cos(2\theta) 阶散光像差
Z_4^4 \sqrt{10}\rho^4\cos(4\theta) 水平四叶草像差

已经选择归一化项,以便

\int_0^{2\pi}\int_0^1 Z_n^m(\rho,\theta)Z_p^q(\rho,\theta)\rho\textrm{d}\rho\textrm{d}\theta=\pi\delta_{n,p}\delta_{m,q}

 

如果两个下标相等,则 Kronecker deltaδ 等于 1,否则为 0。

与我们之前的术语保持一致,泽尼克多项式在权函数 w(\rho,\theta)=\rho 下的单位圆上是正交的。在笛卡尔坐标和圆柱极坐标 \textrm{d}\Omega = \rho\textrm{d}\rho\textrm{d}\theta 之间转换时,权函数 ρ 可以很方便的与雅可比行列式完全匹配。

泽尼克多项式广泛用于光学界。如果你曾经看过眼科医生,那你可能听说过诸如“散光”(Z_2^{\pm 2})) 或“慧差”(Z_3^{\pm 1})之类的术语,而“近视”和“老花眼”只是“散焦”(Z_2^0)一词的不同表述。

一个图形显示了 21 个图的金字塔,将 Zernike 多项式可视化到五阶
高达五阶的 Zernike 多项式图。

用泽尼克多项式表示变形镜

最后,让我们用所学的知识来做一个例题。示例中的几何体是一个侧面和底部固定的平面柱面镜,而上表面可以自由变形。透镜被均匀加热,产生热应力,并导致镜面的上表面膨胀。

代表未变形镜子的短圆柱几何体
未变形镜()和变形镜(右)的几何形状。

我们希望计算最适合拟合变形表面镜位移场的泽尼克多项式系数。如以下屏幕截图所示,为了完整起见,我们将包括所有四阶泽尼克多项式,尽管只有径向指数为 0 的项(平移 Z_0^0, 散焦 Z_2^0, 球差 Z_4^0) 由于位移场的对称性而发挥了重要作用。为了对位移与每个泽尼克多项式的乘积进行积分,在变形表面上定义了积分 耦合。然后使用一些变量 节点来定义每个泽尼克多项式与解向量的乘积的积分,这决定了相应泽尼克系数的值。

最后,我们绘制出位移场w在从镜面中心向外径向延伸的切线上的面外分量,并将其与泽尼克多项式的线性组合进行了比较。下图显示了平移、散焦和球差的单独组合,并将它们的线性组合与位移场进行了比较。解与泽尼克多项式拟合的一致性看起来相当好,并且可以通过引入高阶项(如Z60)进一步改进。

COMSOL Multiphysics 中比较位移场的线图,蓝色显示,泽尼克多项式拟合紫色,活塞为绿色,散焦为红色,球差为青色
位移场(蓝色)与泽尼克多项式拟合(紫色)比较。分别显示了平移(绿色),散焦(红色)和球差(青色)

结论

曲线拟合是一种消除实验数据中的统计噪声的工具,它也可以应用到模型求解本身,效果很好。将解数据拟合到一组预定义函数,特别是正交函数中,提供了一种快速方便的数据压缩方法。

下一步

通过点击下面的按钮来尝试这篇博客文章中讨论的模型,将链接到COMSOL案例下载条目和可下载的mph文件。

参考文献

  1. ISO 24157:2008: Ophthalmic optics and instruments — Reporting aberrations of the human eye, International Organization for Standardization, Geneva, Switzerland. Amendment 1, ibid., 2019.
  2. ANSI Z80.28-2017: American National Standard for Ophthalmics — Methods of Reporting Optical Aberrations of Eyes. American National Standards Institute, Alexandria, VA.

 

博客分类


评论 (1)

正在加载...
tx 孙
tx 孙
2021-10-12

请问如何在绘图组中导出空间间距一致的二维/三维绘图数据呢

浏览 COMSOL 博客