Fenglielie

关于 Euler 方程组的一些常用的基本知识。 控制方程与状态方程 一维的可压缩流体满足 Euler 方程组，主要由三个控制方程组成 \[ \begin{pmatrix} \rho \\ \rho u \\ E \end{pmatrix}_t + \begin{pmatrix} \rho u \\ \rho u^2 + p \\ u(E + p) \end{pmatrix}_x = 0 \] 分别对应质量守恒、动量守恒和能量守恒。涉及的记号含义为： 密度（density） \(\rho=\rho(x,t)\)； 速度（velocity）\(u=u(x,t)\)； 压强（pressure）\(p=p(x,t)\)； 单位体积的总能量（total energy）\(E=E(x,t)\)。 这里还缺少一个方程与控制方程组一起组成封闭的方程组。 可压缩流体的单位体积总能量\(E\)包括动能和内能两部分 \[ E = \frac12 \rho u^2 + \rho e \] 其中\(e\)是比内能（specific internal energy），由流体的压强和密度决定：\(e =...

概述 在学习 C/C++ 的一开始，我们就接触到基本的输入输出，例如printf/scanf或者cout/cin， 但是在涉及用户输入时总是存在可靠性问题：我们无法确保用户输入的信息是合法的，比如要求输入一个数字，但是用户提供的是一串字符串等等，或者要求数字满足一定范围等，我们只能每次手动在外层加一个循环和判断条件，非常繁琐。 在本文中实现了一个简单的获取安全输入的封装SafeInput，使得在编程中尽可能摆脱这类涉及用户输入参数的处理细节。 使用示例 指定获取一个int类型的值，要求大于 0 且小于 5，如果输入不满足要求，则会提示重新输入 12int s = SafeInput().get<int>("Please input a int in (0,5): ", [](int p) { return 0 < p && p < 5;...

本文包括常见的高阶的有限差分公式，以及使用 Mathematica 自动计算数值微分公式的代码。 原理就是基于待定系数法求解得到如下形式的有限差分公式 \[ f^{(\alpha)}(x) = \frac{1}{h^\alpha}\left( c_{-m} f(x-m h) + c_{-m+1} f(x-(m-1)h) +\cdots + c_n f(x+n h) \right) + O(h^\beta) \] 其中的\(m,n \ge 0\) 代表模板宽度，要求\(m+n \ge \alpha\)。 Mathematica 源代码 1234567891011121314FDFormula[m_, n_, a_, {f_, x_, h_}] := Module[{mat, b, coeff, expr, i, j}, mat = Table[If[j == 0, 1, i^j], {j, 0, m + n}, {i, -m, n}]; b = Table[If[i == a, a!, 0],...

任意拉格朗日法是一种用于处理动边界问题的方法，尤其在流体问题中得到应用，以这个方法的学习作为这一系列笔记的结束。 参考文献如下 Donea, Jean, Antonio Huerta, and A Rodrıguez-Ferran. “Chapter 14 Arbitrary Lagrangian–Eulerian Methods,” n.d. Fehn, Niklas, Johannes Heinz, Wolfgang A. Wall, and Martin Kronbichler. “High-Order Arbitrary Lagrangian-Eulerian Discontinuous Galerkin Methods for the Incompressible Navier-Stokes Equations.” Journal of Computational Physics 430 (April 2021): 110040. https://doi.org/10.1016/j.jcp.2020.110040. 1. 不可压 NS 方程 在区域...

关于 Gauss 数值积分公式的推导以及相关性质的笔记，注意并不是与正态分布相关的高斯积分 \(\int e^{-x^2}\,dx\)，在 wiki 上 Gauss 数值积分公式又称为 Gauss 求积公式。 在本文中使用的下标范围为 \(0-n\) 而非 \(1-n\)，这样与多项式次数更加对应，但是与某些资料上的公式不一致。 Gauss 积分推导 最优代数精度 对于标准积分区间 \([-1,1]\) 的数值积分公式，节点和系数的选取有 2n+2 个自由度，至多可以达到 2n+1 阶的代数精度。 \[ I(f) = \int_{-1}^1 f(x)\,dx \approx I_n(f) = \sum_{i=0}^n A_i f(x_i) \] 断言：无论节点和系数如何选取，上述形式的数值积分公式不可能具有 2n+2 阶的代数精度。 证明断言：构造 2n+2 次多项式 \(p(x)\) \[ p(x) = \Pi_{i=0}^n (x-x_i)^2 \ge 0 \] 显然有 \[ \begin{aligned} 0 < \int_{-1}^1 p(x)\,dx =...

主要内容是最优化中的牛顿法和拟牛顿法，还包括最简单的最速下降法，都是用于求解无约束最优化问题的梯度类最经典算法。 对于无约束最优化问题 \[ \min_{x \in R^n} f(x) \] 通用的算法流程如下： 初始化，选取起点 \(x^0\)，记 \(k=0\) 在当前位置 \(x^k\)，获取搜索方向 \(d^k\)，通常是基于当前位置的一阶导和二阶导信息 在给定方向上进行一维搜索，\(a^k = \text{argmin}_{a \ge 0} f(x + a d^k)\)，获取步长因子 \(a^k\) 更新 \(x^{k+1} = x^k + a^k d^k\)，\(k=k+1\)，回到第二步 关于一维搜索的细节，这里不做讨论，这篇笔记只关注搜索方向的确定。 最速下降法 每一次的搜索相当于对 \(f\) 进行了一阶的泰勒展开。 \[ \begin{aligned} f(x) &= f(x^k) + \nabla f(x^k)^T (x-x^k) + O(|x-x^k|^2) \\ \end{aligned} \] 这里记...

本文只考虑最简单的一维问题，包括求根和求最小值点两个子问题： 非线性方程求根 \(f(x)=0\) 非线性方程求最小值点 \(x_0 = \text{argmin} f(x)\) 精确一维搜索 非精确一维搜索 由于线性问题是平凡的，因此默认\(f(x)\)是非线性的。 非线性方程求根 对分法 假设\(f(x)\)在有限区间\([a,b]\)满足：\(f(a)f(b)<0\) 异号，那么\(f(x)\)在\((a,b)\)至少有一个零点 \(x^*\)。 算法： 起始的计算区间 \([x_l,x_r]=[a,b]\)，此时满足两端异号 \(f(x_l)f(x_r)<0\) 。 选取中点 \(x_m = \frac12(x_l+x_r)\) ，计算\(f(x_m)\)。 如果中点 \(|f(x_m)| \le \varepsilon\)，或者区间长度\(|x_r-x_l|\le \varepsilon\) 已经足够小，则认为中点已经足够精确，已经找到零点，求解完成。 否则左右两个区间...

古典迭代法可以用于求解线性方程组，特点是格式简单，收敛性与迭代初值无关，与迭代矩阵的谱半径有关。 古典迭代法实质上是基于矩阵的近似逆设计的不动点迭代法，与此不同的，共轭梯度法等则是最优化在矩阵求解问题上的具体体现。 1. 迭代格式与收敛性 求解非奇异的线性方程组 \[ A x = b \] 将方阵\(A\)分解为三个部分 \[ A = D-L-U \] \(D\)为对角阵 \(L\)为下三角阵，对角线全零 \(U\)为上三角阵，对角线全零 以下的几类迭代法都有如下形式（实质即一种不动点迭代法） \[ x^{(k+1)} = B x^{(k)} + g \] 最终希望收敛到 \(Ax=b\) 的精确解 \(x\)（在向量范数意义下） \[ \lim_{k\to \infty} \Vert x^{(k)} - x \Vert = 0 \] 此时\(x\)满足如下等价形式 \[ x = B x + g \] 收敛的充要条件与 \(B\) 的谱半径有关： \[ \rho(B) < 1 \] 收敛的充分条件有很多，例如 \(B\) 的某个矩阵范数小于 1...

一维 Burgers 方程如下 \[ \left\{ \begin{aligned} u_t + \left(\frac{u^2}2\right)_x &= 0\\ u(x,0) &= u_0(x) \end{aligned} \right. \] 其中定义域 \(x \in \mathbb{R}\)，通常假定初值有周期性，精确解此时也具有周期性。 特征线方程 从\((\xi,0)\)发出的曲线记作 \(l_\xi\)，曲线参数为 \(s\)，方程为 \(x = x(s)\)，其中 \(x(0)=\xi\)。沿着曲线 \(l_\xi\)，\(u\) 的变化情况为 \[ \begin{aligned} \frac{d}{ds}u(x(s),s) = u_x(x(s),x) x'(s) + u_t(x(s),s)\\ = u_x(x(s),s)\, (x'(s) - u(x(s),s)) \end{aligned} \] 若曲线 \(l_\xi:x=x(s)\) 取为如下...

主要参考了Google test(gtest)和知乎上的一篇文章qtest: 一个单元测试库的从头实现以及作者提供的代码，尤其是宏的部分。一直不喜欢也没有学明白宏的各种用法，但是实现这种风格的测试框架也绕不开宏。 在其基础上进行了整理和重构，并且扩展和完善了一些细节的功能。 本文作为 C++ 学习中的小练习，若有疏漏，欢迎指正。 MTest 介绍 MTest 是一个 header-only 的简易框架，只包含两个头文件： mtest.hpp 负责 MTest 和 MTest::MTestMessage 两个类的实现，不含有任何的宏。 mtest_macro.hpp 负责对外提供相应的宏，可以在编译时使用-DUNUSE_MTEST关闭所有的宏，避免与 gtest 产生冲突。 测试文件需要 include 这两个文件，include顺序无所谓，也可以直接合并成一个文件，但是我个人的喜好是对宏敬而远之，因此单独仍在一个头文件中。主要功能实现在mtest.hpp，其中不含有任何的宏。 为什么重复造轮子？ 我希望将 MTest 作为 gtest 的简易替代，尤其是它具有...