智障。 - 数值计算
type: channel
link: http://www.dmcxe.com/index.php/tag/数值计算
items: 数值计算三——微分方程组的求解与Modelica
items: 传热学二：利用稀疏矩阵快速求解
items: 数值计算一——插值、拟合与线性方程组

数值计算三——微分方程组的求解与Modelica
type: item
date: 2023-04-22T10:40:03+00:00
description: 本次将通过改进RK4（四阶龙格-库塔法）法构建ODE求解器以获得对于微分方程组的求解，并通过相关大作业的实力内容进行计算验证。通过此方程组的建立，对于Modelica运行时的一些event耗时问题也产生了一些思考与理解。此外，代码全部上传至https://github.com/DMCXE/Numerical_ComputationODE.py Homework.moRK4方法的回顾常用的ODE求解方法为龙格-库塔方法。即增加欧拉法中的节点数提高精度。RK方法可以用公式表示为： y n + 1 = y n + h φ ( x n , y n , h ) , 其中 φ ( x n , y n , h ) = ∑ i = 1 r c i K i K 1 = f ( x n , y n ) K i = f ( x n + λ i h , y n + h ∑ j = 1 i − 1 μ i j K j ) , i = 2 , ⋯ , r 当对r取不同的取值，则对应了不同的方法。四阶RK4方法具备较高的精度，常用于数值计算中。通过数学计算，可以得到RK4的一般表达式为： { y n + 1 = y n + h 6 ( K 1 + 2 K 2 + 2 K 3 + K 4 ) , K 1 = f ( x n , y n ) , K 2 = f ( x n + h 2 , y n + h 2 K 1 ) , K 3 = f ( x n + h 2 , y n + h 2 K 2 ) , K 4 = f ( x n + h , y n + h K 3 ) . 但通常微分方程是通过微分方程组的形式出现，比如说传染病模型SEIR。而在《数值计算二》中构建的ODE求解器只能够求解单独的微分方程，局限性较高。开发通用的微分方程组求解器是必须需要解决的。RK4方法对方程组的适应改进在这里，RK4是一种前向方法。因此当函数化为方程组时，只需要将方程指定为向量格式、并在自变量统一时，将因变量写成向量格式。若f(x,y)表示的形式是f(t,x)，即写为f(t,x1,x2,x3,x4)。RK4的数学形式是支持向量化的，因为按照上述方法处理后，各个变量的相互关系能够通过向量组缩影体现出来。在求解器中，求解过程完全一致，但与求解单方程不同的是，我们需要指定方程组在的方程个数N，并通过该数构建用于储存结果的Nx1矩阵。我们把主要的处理放在了对于函数组表达的预处理上。import numpy as np class RK4_for_equations: def __init__(self, F,n,min,max,totstep,begin): self.F = F self.n = n #指定方程组的个数 self.min = min self.max = max self.step = (self.max-self.min)/totstep self.begin = begin self.totstep = totstep def slover(self): f = self.F x = 0 y = self.begin yn = np.zeros((self.n,1)) yn[:,0] = y step = self.step flag = 1 for i in range(1,self.totstep): K1 = f(x,y) #print("K1:",K1) K2 = f(x+0.5*step,y+0.5*step*K1) #print("K2:",K2) K3 = f(x+0.5*step,y+0.5*step*K2) #print("K3:",K3) K4 = f(x+step,y+step*K3) #print("K4:",K4) y = y + (step/6)*(K1+2*K2+2*K3+K4) yn = np.append(yn,y.reshape(self.n,1),axis=1) flag = flag + 1 x = i * step return yn对于传入函数的预处理简单的来说，我们需要把传入函数写为F(x,y)=[f1(x,y,...,fn(x,y)]格式。以下通过一个实际案例进行描述：已知蒸汽发生器水位控制系统的微分方程如下： x ̇ 1 = G 1 ( q e − q v ) x ̇ 2 = − 1 τ 0 x 2 + G 2 τ 0 q v x ̇ 3 = 1 τ 2 x 4 − G 3 β τ q e x ̇ 4 = − x 3 − 2 ξ τ x 4 + ( 2 ξ β − 1 ) G 3 q e y ( t ) = x 1 + x 2 + x 3 x 1 ( 0 ) = x 2 ( 0 ) = x 3 ( 0 ) = x …

传热学二：利用稀疏矩阵快速求解
type: item
date: 2023-01-19T11:40:46+00:00
description: 本次介绍直接通过求解线性方程得到二维平板的稳态传热。当然，醉翁之意不在酒，拉普拉斯方程代表的一众物理问题都可以迎刃而解。本内容上传至学院共享仓库https://github.com/euaurora/HappyShares-CNST-HEU/tree/main/传热学/平板传热的矩阵法%20-by%20DMCXE以及个人博客。推荐观看pdf版本。1 回顾严格来说是前年，介绍了如何通过Gauss–Seidel迭代方法求解二维平板传热方程。以Gauss–Seidel、Jacobi等为代表的一种方法均通过迭代进行矩阵的求解。由于其收敛曲线平滑，因而对于满足此类方法的矩阵，复杂与否都可以得到期望的结果。但是缺点显而易见，在上一篇文章中，100X100的网格计算了90s，这显然是不可接受的。对于传热方程、或者对于形式上类似的方程， 我们有更好的方法进行求解。2 矩阵运算快速求解2.1 二维传热方程的矩阵运算无内热源，稳态的二维传热方程为: ∇ 2 ϕ = 0 空间差分可以得到： ϕ i − 1 , j − 2 ϕ i , j + ϕ i + 1 , j Δ 2 x + ϕ i , j − 1 − 2 ϕ i , j + ϕ i , j + 1 Δ 2 y = 0 考虑均匀的网格划分$\Delta x= \Delta y$, 上述方程可以化简为: ϕ i − 1 , j + ϕ i + 1 , j + ϕ i , j − 1 + ϕ i , j + 1 − 4 ϕ i , j = 0 根据《传热学》，这一步就可以进行迭代表达式的拆分了。直观来看，似乎并没有什么好方法去求解这个复杂的线性方程组。2.1.1 在一维情况下思考如果我们退一步，考虑一维的传热方程 ϕ i − 1 − 2 ϕ i + ϕ i + 1 = 0 我们遍历每一种可能情况： ϕ 1 − 2 ϕ 2 + ϕ 3 = 0 ϕ 2 − 2 ϕ 3 + ϕ 4 = 0 ϕ 3 − 2 ϕ 4 + ϕ 45 = 0 注意到可以写为一个系数矩阵A_{NxN}与T_{Nx1}的乘积，且系数矩阵只有三个对角！（不考虑周期边界条件） A = [ − 2 1 0 ⋯ 1 1 − 2 1 ⋯ 0 ⋮ ⋮ ⋱ ⋯ ⋮ 0 ⋯ 1 − 2 1 1 ⋯ 0 1 − 2 ] N × N T = [ ϕ 1 ϕ 2 ⋮ ϕ N ] Z = [ 0 0 ⋮ 0 ] 也就是说，传热方程可以写成一个线性方程组AT=Z。求解矩阵的线性方程组，自然而然的就想到几个方法：求逆。T=A^{-1}Z。分解。A=LU,LUT=Z,b=UT,Lb=Z也就是说，这个方程好不好解，除了与本身的规模大小有关外，还与系数矩阵A的性质有关。我们重新观察A，发现他是：稀疏的。只有几个对角。——> 稀疏矩阵的逆往往是稠密的，稠密矩阵进行运算复杂度是N^3对称的。A=A^{T}在有边界值的情况下，是非奇异的。A的性质非常好，也就是说可以通过特殊的矩阵分解方法进行快速矩阵运算求解。举个例子，加入对于一维杆，两端温度已知，上述矩阵为： A = [ 1 0 0 ⋯ 0 1 − 2 1 ⋯ 0 ⋮ ⋮ ⋱ ⋯ ⋮ 0 ⋯ 1 − 2 1 0 ⋯ 0 0 1 ] N × N T = [ ϕ 1 ϕ 2 ⋮ ϕ N ] Z = [ T 0 0 ⋮ T N ] A构成了三对角矩阵，求解三对角矩阵的线性方程组可以使用追赶法：对于方程组 ( b 1 c 1 a 2 b 2 c 2 ⋱ ⋱ ⋱ a n − 1 b n − 1 c n − 1 a n b n ) ( x 1 x 2 ⋮ x n − 1 x n ) = ( f 1 f 2 ⋮ f n − 1 f n ) 利用1.计算{beta_i}的递推公式： β 1 = c 1 / b 1 β i = c i / ( b i − a i β i − 1 ) , i = 2 , 3 , ⋯ , n − 1 2.求解Ly=f y 1 = f 1 / b 1 y i = ( f i − a i y i − 1 ) / ( b i − a i β i − 1 ) , i = 2 , 3 , ⋯ , n 3.求解Ux = y x n = y n x i = y i − β i x i + 1 ， i = n − 1 , n − 2 , ⋯ , 2 , 1 即可实现快速计算，即一种特殊LU分解。2.1.2 在二维中回到 ϕ i − 1 , j + ϕ i + 1 , j + ϕ i , j − 1 + ϕ i , j + 1 − 4 ϕ i , j = 0 如果假设我们存在一个系数矩阵A_{NxN}$与描述平面温度分布的矩阵T_{NxN}，这样的稠密矩阵相乘是非常繁琐的。自然的，类比一维空间差分的处理，自然的想到该方程对于…

数值计算一——插值、拟合与线性方程组
type: item
date: 2022-10-11T13:03:00+00:00
description: 这次，我们将通过插值、拟合与线性方程组的数学描述导出数值方法，并用其解决研究生课程《数值计算》大作业中的部分问题。本文涉及的全部代码及其用法可前往Github仓库https://github.com/DMCXE/Numerical_Computation 查看。目录如下： 我强烈建议您使用Github阅读本文章，因为我写完才想起来目前无法适配Latex公式！！！全文共1w字，至少3000都是公式！！！仓库地址：https://github.com/DMCXE/Numerical_Computation我在Github的公式表达上遇到了繁琐的问题，在彻底修复前可以参阅pdf文件！1 插值法我们为什么要插值？在一个核裂变反应堆内，为了精细化建模我们需要通过试验确定不同核素和不同能量中子发生相互反应时的反应截面。最好的情况是，有一套完整的作用理论，我们只需要测几个普遍数据作为验证即可。但我们不仅尚未构建出截面模型，甚至依托实验手段无法精准的得到不同能量得到中子束以测定不同能量下的反应截面。更致命的是，由于多普勒（Doppler）效应，目标核素自身热运动大小改变了中子的相对能量，使得截面产生多普勒展宽（Doppler-Broadening）。运用插值法，我们希望解决这样的问题：通过有限的能量-截面数据推出任意能量下的截面数据。确保推导产生的数值是精确的，符合客观规律的。产生如多普勒展宽效益下，局部插值同样能发生改变，且保持其余区域的稳定性快速计算，节省计算资源，提高分析效率我们希望根据给定的函数表做一个既能反应截面函数$f(x)$的特性，又便于计算的简单函数$P(x)$，用$P(x)$近似$f(x)$。而这便同样是插值法的目的。1.1 插值方法的定义设函数$y = f(x)$在区间$[a,b]$上有定义，且已知$a \leqslant x_0<x_1<\cdots<x_n \leqslant b$上的值$y_0,y_1,...,y_n$，若存在一简单函数$P(x)$，使：$$P(x_i)=y_i, i = 0,1,...,n$$成立，就称$P(x)$为$f(x)$的插值函数，点$x_0, x_1, \cdots, x_n$称为插值节点，包含插值节点的区间$[a,b]$称为插值区间，求插值产生$P(x)$的方法称为插值法。依据插值函数的不同，分为：插值多项式：$P(x)$是次数不超过n的代数多项式分段插值：$P(x)$为分段多项式三角插值：$P(x)$为三角多项式1.2 多项式插值定理：设在区间$[a,b]$上给定的$n+1$个点$a \leqslant x_0<x_1<\cdots<x_n \leqslant b$上的函数值$y_i = f(x_i)(i=0,1,\cdots,n)$，求次数不超过$n$次的多项式，使得$$ P(x_i)=y_i, i=0,1,\cdots,n $$满足此条件的多项式是唯一的。1.2.1拉格朗日（Lagrange）插值多项式1.2.1.1 数学定义定义：若n次多项式$l_j(x)(j=0,1,\cdots,n)$在$n+1$个节点$ x_0<x_1<\cdots<x_n$上满足条件：$$ l_j\left(x_k\right)=\left\{\begin{array}{ll} 1, & k=j, \\ 0, & k \neq j, \end{array} \quad j, k=0,1, \cdots, n,\right. $$就称这$n+1$个$n$次多项式$l_0(x), l_1(x), \cdots, l_n(x)$为节点$ x_0<x_1<\cdots<x_n$上的n次插值基函数。可以表示为：$$ L_n(x)=\sum_{k=0}^n y_k l_k(x) $$其中，$l_k(x)$表示为：$$ l_k(x)=\frac{\left(x-x_0\right) \cdots\left(x-x_{k-1}\right)\left(x-x_{k+1}\right) \cdots\left(x-x_n\right)}{\left(x_k-x_0\right) \cdots\left(x_k-x_{k-1}\right)\left(x_k-x_{k+1}\right) \cdots\left(x_k-x_n\right)}, \quad k=0,1, \cdots, n . $$满足上式的多项式可称为拉格朗日(Lagrange)插值多项式。不失简便性，可写成如下形式：$$ L_n(x)=\sum_{i=1}^n y_i\left(\prod_{j \neq i}^{1 \leq j \leq n} \frac{\left(x-x_j\right)}{\left(x_i-x_j\right)}\right)…