https://www.youtube.com/playlist?list=PLwXscOMmh_DJaSIE_RtDO2FDyey5wAkox
该课程主要包含差分法、有限体积法、有限元素法和积分边界元素法四种求解偏微分方程的数值方法。还将重点讨论如何将偏微分方程 discretization 成系统方程,然后使用直接求解方法解决离散后的方程。
偏微分方程广泛用于天气预测、航天飞船设计以及建筑消防安全规定等领域。数值模拟可以替代实验,从而节省成本和时间。
学生提出希望学习不同方法的限制条件,如何处理非线性问题,以及如何根据问题选择适当方法。重要的是掌握编写代码求解方程的能力。
将按顺序学习差分法、有限体积法、有限元素法和积分边界元素法。同时重点授课直接求解离散后的系统方程的方法。帮助学生养成解决实际问题的能力。
该课程没有考试,但将重点授课如何解决实际问题。为此,课程将设置大量项目,要求学生编写代码并调试以解决给定问题。
第一个个人项目会评分,但不计入总成绩。之后的项目将以团队形式进行,首个团队项目成绩将计入总成绩。
项目中禁止在团队或个人之间共享任何文档或代码。允许讨论问题,但结果不得记录在文件或代码中。
本课程将主要使用MATLAB进行编程。也允许使用Python,但当前学生MATLAB水平较高,所以目前采用MATLAB。
需要在规定时间内完成第一个个人项目,以便给予合适的团队组合。团队项目结果将计入总成绩。
本课程使用一个泮疸方程作为实例来介绍数值解法:
ut + uux = κuxx
该方程在工程中广泛应用,可以描述温度的传播、污染物浓度的蔓延以及油气田中的水浸现象等。
该方程是线性方程,因为如果u1和u2均为解,则u1+u2也是解。
给出初值条件u(x,0)=u0(x)和边界条件u(0,t)=uL(t), u(L,t)=uR(t)。
使用MATLAB模拟方程数值解,验证κ、u各参数组合下解行为。
理解方程的物理意义有助于验证数值算法是否正确。这也是选择该方程的原因。
当对流项u等于0时,泮疸方程会退化为热方程:
ut = κuxx
热方程描述固体或流体的温度随时间和空间的变化,可以模拟初始温度分布随热传导而逐渐均匀的过程。
热方程属于泮疸型方程,其特征是初始条件的变化将在有限时间内影响整个计算域。
用有限差分法和有限体积法对连续函数进行 discreization,以在计算机中对函数进行离散表示。
高阶方法是指当网格细分程度加倍时,解的误差下降速度大于线性。一般研究第一阶和第二阶方法,也会学习第三阶以上的高阶方法。
解释有限差分法和有限体积法的区别,以及高阶方法与 discreization 方法的关系。
有限差分方法将连续函数近似表示为计算点的值,仅保留函数在空间格点上的离散表示。将热方程问题域离散为均匀间隔的格点。
首先在空间上离散,再将函数表示为向量对时间的函数。将偏微分方程转换为一组常微分方程,获得时间上的离散表示。
采用Forward Euler方法解常微分方程组,利用已知初始条件按时间步长递推求解。
利用泰勒展开近似表达邻点函数值,将空间导数近似表示为函数值差 quotient,取得有限差分式。其误差随格点间隔缩小而缩小。
有限差分法近似误差的顺序决定了扩网后误差下降的速率。一般采用二阶或更高顺序方法取得更好结果。
Matlab提供影响微分方程解的函数,可以直接求解离散化后的微分方程组,获得网络格点上的离散解。
有限差分求解时需要考虑边界条件。如果已知边界函数值,相应边界点的解直接给定;否则边界点的时间微分与相邻内点关系式会改变。
将有限差分求解表达为ODES矩阵形式du/dt=Au+B,其中A为三对角矩阵,B为考虑边界条件后修正项。
利用Forward Euler和Backward Euler对时间微分进行离散,将问题转化为求解线性方程组。Forward Euler为显式法稳定性差,Backward Euler稳定性好但需要求解方程组。
Forward Euler和Backward Euler均为一阶时间近似,其时间步长对数值误差的影响较小。PDE求解误差与时间近似顺序和空间近似顺序同时影响。
如Trapezoidal法为二阶时间近似,Runge-Kutta法顺序更高。MATLAB中ode45采用自适应四阶或五阶Runge-Kutta法获得更高精度。
教师发布了Project 0的要求,要求学生实现已讨论过的热方程的有限差分数值求解。具体来说:
自行推导有限差分方案,包括边界条件;
写一个MATLAB函数heat1D.m,输入步长Δt,边界条件和初值函数,输出域中间位置的解。
要求提交PDF文件,不包含源代码。源代码单独文件。
文件清晰扫描后提交,以便评估同学们的编程水平。
这次作业不作为期末成绩考核,主要目的是帮助教师了解课堂人数,以便搭配合作伙伴。
可以手写或键入结果;
扫描成PDF文件时,放置于光线不直接的位置;
仅听课的学生,可以与其他听课生搭档完成。
上次课讨论了一维热方程的有限差分解法,将导数近似为迹差商,将方程离散为时空网格点上的线性方程组。
类似地,可以将二阶偏微分方程如Δu=f离散为网格点上的线性系统Au=F。
与熱方程相比,波义森方程不含时间,直接求解线性系统取得steady state解。
两方程离散后皆为时间的向量微分等于矩阵乘以解向量,或线性系统的形式。
但熱方程用ode45求解,波义森方程直接用MATLAB \或inv(A)*F求解。
下一步将讨论FD近似的误差性质,关键是稳定性条件。也将介绍如何自行构造FD方案。
将提到基于同样原理可扩展到二维情况,离散为网格点系统即可。
微分方程∂^2u/∂t^2=∂^2u/∂x^2描述一维域内波动传播。其解u(x,t)可以是x-t或x+t的函数,表示波向两方向传播。
撮扰在时空点只影响有限范围内,形成两边向外扩散的锥面。信息传播速度有限,属于超杆方程。
将∂^2u/∂t^2代表为新变量v=∂u/∂t,得到u_t=v,v_t=∂^2u/∂x^2成为u,v随时间的耦合微分方程。
使用同样的差分格式近似Spatial微分,将方程离散为动量U,V在网格点间的耦合变化,用矩阵表达。采用OD45等数值求解。
超杆方程的通性是信息传播具有有限速度,形成锥状域影响范围。其形状可能不同,但求解方法相同。
截断误差定义为数值差分近似与实际微分值之间的差别。它反映了数值近似的精确程度。
通过泰勒展开计算截断误差。将UI+, UI-, UI展开至足够低阶,保留第一个非零项。
将各项求和,除以Δx^2。除第一阶导数外,所有低阶导数项互相抵消。留下Δx^2/12×第四阶导数项。
截断误差取决于Δx和解函数的第四阶导数。Δx越小,误差越低。振荡函数高阶导数大,使用同样Δx时误差也大。
要使不同函数截断误差一致,需为振荡函数使用更小的Δx。可降低误差或者提高解空间的分辨率。
通过数值近似方式求解导数,并与真实解进行对比,得到截断误差。可用泰勒展开求解真实解,找到预测误差大小的最低阶项。
通过绘制数值解与预测误差之间关系的对数曲线来检验数值实现是否正确。若曲线为直线,且斜率等于误差阶数,则实现正确。
Δx趋于0时,误差近似为Δx的误差阶数次方。每减小Δx一个数量级,误差应减少两个数量级。这可以用于判断数值实现是否达到预期误差阶数。
前面是前验误差分析,需要真实解。后验误差分析通过数值解写出误差上限,不需要真实解,但难度大。给出误差与数值解的关系,通常为不大于关系。
泰勒展开可以寻找截断误差的数学表达式。通过绘制误差对数曲线检验数值实现,判断误差随Δx变化趋势,评估方法有效性。后验误差分析给出通过数值求得的误差上限。
通过泰勒展开,比较近似表达式与真实表达式之间的差异来构造有限差分近似公式。
构造公式需要满足:
只采样点的直接相邻点,求得紧密简洁的表达。
截断误差应达到最低阶,通常为Δx的次方。
设置逼近表达式为点的线性组合,引入未知系数。
通过泰勒展开确定近似表达式与真实表达式之间的关系。
设置表达式中未知项为0,消去低阶项,得到线性方程组。
用矩阵法解方程组求解未知系数,得到最终近似公式。
如 central difference 方案就是此方法得到的二阶有限差分近似。
边界点采用不同方案:1. 边界点设计特殊方案;2. 利用边界条件约束方案。
除中点差分外,前后差分方案也是一阶近似,但在求解时稳定性可能较好。
在解析有限差分方案的误差时,需要区分截断误差与数值解误差的含义:
截断误差为有限差分近似算子与行列式求导算子在解析解上的差异。
数值解误差为数值解与解析解之间的差异。
将行列式方程组记为Au=f,其中A为有限差分近似算子,u为数值解,f为给定函数。
则Au-Af=A(u-u)=数值解误差称为截断误差。
数值解误差与截断误差之间的关系为:
数值解误差=A^-1×截断误差
其中A^-1为A的逆矩阵。
要求数值解具有一致性,需要满足:
截断误差随Δx收敛至零
矩阵A^-1的范数随Δx保持不变
只有在满足上述两个条件时,数值解误差才能与截断误差具有相同的收敛顺序,从而实现一致性。
稳定性是联系一致性与收敛性的关键概念,它要求矩阵A^-1的范数与Δx无关。
保证方案的稳定性是构造高效有限差分方案的重要考虑。
为了分析有限差分方法的稳定性,需要对矩阵A进行谱分解。
A矩阵对应泊松方程的有限差分表示,其形式为:-2/Δx^2,1/Δx^2,1/Δx^2,-2/Δx^2。
A矩阵为对称矩阵,其特征值和特征向量都是实数。
特征向量可以表示为正弦函数形式:Xk=sin(kπi/n)
对Xk应用A矩阵,可以得到特征值λk的表达式:
λk=Δx^2[sin(kπ/n)-2+sin(-kπ/n)]
最小特征值λ1对应k=1,其值不依赖Δx。
矩阵A的反矩阵范数取决于最小特征值λ1。
若λ1保持有界,则A矩阵稳定,有限差分解将保持一致性。
泊松方程的有限差分表示符合上述条件,其稳定性得以证明。这是构造高效数值算法的重要考虑点。
上课总结了一维泊松方程数值解法的稳定性分析:
错误项包括截断误差和解误差。解误差由截断误差及解算误差矩阵A的逆矩阵相关。
A矩阵对应有限差分表示为:-2/Δx^2,1/Δx^2,1/Δx^2,-2/Δx^2。
A矩阵为对称矩阵,其特征值λk可以表示为函数Δx形式。
计算A矩阵逆矩阵范数的上界,取决于特征值λ1。
应用泰勒级数展开泊松方程的有限差分表达式,当Δx趋于0时,A矩阵逆矩阵范数上限为L^2/π^2,与Δx无关。
这意味着解的误差随Δx的减小会以与截断误差同样的速度减小,保证了数值方法的稳定性。
上述分析适用于一维泊松方程,同样方法可扩展至其他偏微分方程。分析核心是求解特征值与特征向量,研究A矩阵逆矩阵的上界。
本节主要讨论一维泊松方程的边界条件:
边界条件主要有Dirichlet条件和Neumann条件两种。
Dirichlet条件就是解函数在边界点的值被固定。
Neumann条件就是解函数在边界点的一阶导数值被固定。
以热传导方程为例,Dirichlet条件对应固定温度的热源,Neumann条件对应固定热流的绝热边或热源边。
对于Neumann边界条件,有限差分矩阵A的第一行会有修改。以一阶精度为例,第一行为[-1,1]/Δx,对应导数近似。
高阶近似需要采用中心差或偏离差格式,以提高截断误差顺序。
也可以考虑混合边界条件,即解函数和其一阶导数的线性组合被固定。
边界条件的选择会影响有限差分方程的解法,正确描述实际问题的边界也很重要。
本节主要介绍二维场景如何使用域映射将问题映射为矩形域,以便于应用二维有限差分方法求解:
场景几何结构复杂时,可以通过映射将物理域映射为计算域(矩形域)。
两个典型映射方式为:锯齿状映射和鱼钩映射。它们可以将飞翼等复杂形状映射为矩形。
由映射得到物理域坐标(x,y)与计算域坐标(ξ,η)之间的转换函数关系。
将 Partial differential equations 从(x,y)坐标转换为(ξ,η)坐标,利用链式法则将导数转换为新坐标下的形式。
转换后方程中会出现映射派生多项式,正确选择映射可以使一些项消失,简化计算。
进而可以在计算域上使用有限差分方法求解转换后的方程,从而获取物理域解。
映射可以精确地对应计算域与物理域的边界,保证解的一致性。正确进行映射是应用二维有限差分的关键。
本节主要回答了域映射方法中的几个问题:
如何将偏微分方程从物理域坐标转换为计算域坐标?利用映射函数和链式法则,将导数转换成新的坐标形式。
转换后方程中会出现映射派生多项式,这可能导致方程从线性变为非线性。但由于映射本身的非线性,转换后的方程仍保持线性。
怎样设计映射函数?无统一方法,需要根据问题域的特点通过试错误进行。
计算域网格是否必须是均匀的?虽然不一定,但对有限差分更友好。也可以在不同方向采用不同分辨率。
映射函数应根据问题情形设计物理域和计算域间的网格分辨率。比如壁面附近需要更细致的网格。
在物理域和计算域坐标下分别对问题进行离散呢?应先从物理域坐标导出计算域坐标下对应的偏微分方程,然后在计算域上进行有限差分处理。
本节介绍二维泊松方程的数值解法:
使用固定Δx和Δy将问题域离散成网格点。函数值只存储在网格点上。
二维拉普拉斯算子可以近似表示为x和y方向二阶导数 suma。
x和y方向二阶导数分别使用五点差分近似。
将二维数组函数值展平为一维向量,以求矩阵形式。
对应五点差分公式,矩阵对角线系数为-(2/Δx^2+2/Δy^2) ,上下左右相邻点系数为1/Δx^2或1/Δy^2。
由于网格结构,矩阵非对称三diagonal,而是五diagonal,其中两对角线相隔n点。
矩阵按块划分,每个块对应一横行,内部为三对角结构,但块与块之间有隔断。
有矩阵表示后,可以使用线性代数方法求解泊松方程,并分析差分式的稳定性等属性。
本段视频使用matlab实现二维泊松方程的有限差分法。
设置网格点个数n=8,m=10;步长Δx=0.1,Δy=0.01。
定义矩阵a,大小为n×m,即80×80。初值为0。
填充对角线,a(i,i)设为-(2/Δx^2+2/Δy^2)。
用双层循环填充上下相邻对角线,a(i+1,j)=1/Δx^2,a(i-1,j)=1/Δx^2。
填充与对角线间隔nentries的对角线,a(i+j×n,j×n)=1/Δy^2,a(i+j×n,i+j×n-n)=1/Δx^2。
打印矩阵后可以看出结构为块三对角 form,每个块大小为8×8,之间存在间隔。
讲师更新了本门课程的作业要求:所有文件需放在根目录,不收取硬拷贝作业。
上节课介绍矩阵A近似离散化二维拉普拉斯算子,将偏微分方程转化为矩阵方程Au=f。
极值型方程直接使用MATLAB解Ax=b,或迭代法。
放散型方程将其视为一阶常微分方程组,使用OD45或隐含法解udt=Au。
波动型方程先将其转化为一阶方程,再使用OD45或隐含法解矩阵方程组du/dt=Auv。
隐含法每步都需要解线性方程,将讨论迭代法解大型矩阵方程。
课后讨论将二阶方程转为一阶的原因是:许多一阶ODE软件库,比如MATLAB的OD45,没有二阶ODE对应的好的解法。
常见的直接法(如LU分解法)因存储成本高而适用于小规模问题。
迭代法如Jacobi方法只作用于矩阵中非零元,不依赖矩阵规模,适用于大规模问题。
Jacobi方法将系数矩阵A分解为三个部分:下三角矩阵L,对角矩阵D,上三角矩阵U。
将线性方程组Au=b分解为D^{-1}(b-Lu-Uu),迭代计算解u的值,直到收敛。
首先设定初值u^0,然后计算u^{i+1}=D^{-1}(b-Lu^i-Uu^i),反复进行,u^i逐渐接近真实解u。
方法是否收敛取决于系数矩阵A的性质。通常初值选为零向量,也可选取上一时刻解。
方法收敛意味着u^i和u^{i+1}差异越来越小,最终等于真实解u,满足线性方程组Au=b。
迭代法优点是只针对非零元计算,不需存取整个矩阵,较适用于大规模问题。
Jacobi迭代法通过对角化系数矩阵A,将线性方程组解解为对角矩阵D的解与其他矩阵L和U的加法组合。
对1D梯度场方程,迭代求解方程采用第一类边值条件。系数矩阵A包含主对角线和次对角线,只与相邻节点相关联。
Jacobi迭代公式为:新值=D逆(右端-L旧值-U*旧值)。其中D逆表示1/对角元,L和U代表次对角线上元素。
对2D梯度场方程,系数矩阵A包含4个方向的相邻节点。Jacobi迭代公式同样运用近邻节点信息求解迭代值。
收敛判断标准是residual,即B-A*u趋近于0。迭代终止条件可以参考residual大小与解的误差范围。
Jacobi迭代法的优点是计算简单,只使用非零元,不依赖矩阵规模。但当对角元为0时,无法直接使用Jacobi法,需通过置换变换修正矩阵。
使用Lenna图像作为2D梯度场方程的精确解,采用该图像的红色分量作为解u。
设置单元格边长Δx=1/129,计算右端项B时,将系数矩阵A的操作转换为邻域内线性组合的离散拉普拉斯算子形式。
定义初值采用全0数组,迭代计算unext。Jacobi迭代公式为unext=1/(Δx^24)(B-u上下左右邻域之和)。
设置1000次迭代,每次迭代将unext赋值给u,更新迭代值。
显示迭代结果可以观察解从零初值逐步还原原图轮廓,色域变化趋于精确解。
Jacobi迭代法可用于数值求解各种等式类型,只需根据求解问题设置适当的初始值、迭代公式和终止条件。收敛速度取决于问题特征。
Jacobi迭代中的边界条件被设置为0,使得迭代域内点与边界点的收敛速度不同。
若边界条件为其他函数值,可在每次迭代结束后更新迭代域边界点值。
若边界条件为法向导数,可根据离敏边界一点的解将边界点更新。
二维数组在需要矩阵操作时需展平为一维数组进行。迭代过程中使用邻域差分可避免此步骤。
解在中间位置误差较大是Jacobi迭代的通性,取决于迭代矩阵的最慢特征向量。
分割迭代域独立求解后,内部点误差形式不变,只依赖域形状,与具体图像内容无关。
迭代次数增加后,解较准确,内外点误差差异不明显。
Jacobi迭代求解误差与数值方案收敛误差定义不同。前者是迭代解与精确解的差,后者是数值解与真实解的差。
将Jacobi迭代方程与原方程相减,可以得到迭代误差的迭代关系表达为迭代矩阵R乘误差向量。
若设初值误差为零向量,则第K次迭代误差可以表示为迭代矩阵R的K次幂乘初值误差。
迭代矩阵R的特征值需满足范数小于1,否则含有特征值大于1部分的特征向量将导致误差呈指数散ivergence。
算法长期迭代后,误差向量可以用特征向量基展开,其中对应特征值接近1的分量将占主导地位。
泊松方程Jacobi迭代矩阵特征值分布与解的光滑程度相关,特征值最大的对应特征向量为最光滑解。
Jacobi迭代解求大规模线性方程组,将矩阵分为对角线部分和非对角线部分。
给定真实解为一张彩色飞机图片,利用二阶中心差分方法将图片离散为大小为1024x1024的线性系统。
定义Jacobi迭代函数,输入为初值解和右侧项,输出为下一步迭代解。内实现逐节点更新计算。
对角线项取负值,非对角线项取边节点值之和。每次迭代得到新的解值,逐步求解线性系统。
计算迭代解与真实解误差,多次迭代后误差呈现平滑程度较高的模样,与真实解特征一致。
可视化迭代解与误差,明显看出迭代解逐步接近真实解,误差随迭代收敛减小。
Jacobi迭代函数可重复调用实现多次迭代,求解大规模线性系统。采用该方法可以有效实现图片的数值解析求解。
使用Jacobi迭代求解大小为1024×1024的线性系统,得到的迭代解边界区域较亮,不是理想结果。
检查代码,发现在进行迭代计算时,迭代变量i的取值范围应为2到n,而非原代码中的2到n-1。修改后问题未解决。
重新从零开始,对初始解u赋值全0矩阵,多次迭代仍得到相同结果。
检查右端项b的生成过程,原代码中的迭代变量j的取值范围也应为2到n。修改后b正确生成。
重新对u赋值全0矩阵,正确生成b后,Jacobi迭代可得到理想结果,边界也呈现较暗颜色。
判断是原解u还是B中使用的变量问题,应使用真实解u_exact而非迭代解u。修改后B正确,迭代结果规范。
删除中间变量,防止错误传播。正确设置迭代变量范围和复制变量,即可通过Jacobi迭代实现线性系统求解。
Jacobi迭代可以使求解逼近理想解,得到的迭代解将会非常光滑。这是因为在迭代过程中,偏微分矩阵D与拉普拉斯算子L+U的分割,使得理论解也满足迭代方程。
将迭代误差定义为实际解与迭代解的差异,可以得到误差方程中的迭代矩阵R。R的主特征值决定了迭代次数K后误差形式。
对2D泊松方程,迭代矩阵R对应均值算子。通过特征分析得知,R的特征向量为正弦函数,特征值为函数K1和K2。
当K1=K2=1时,特征值为近1,对应最光滑模式。当K1=K2=N-1时,特征值为近-1,对应最大振荡模式。
Jacobi迭代推导出解趋于光滑,但存在最大振荡模式难 convergent。Multi-grid算法利用解光滑特性提高迭代效率。
实际迭代初解波动较小,难看出最大振荡模式。若初解为随机噪声,后续迭代可以看出此模式。Overall Relaxation可以确保解 remains光滑。
放松Jacobi迭代可以确保迭代解匀称光滑。这通过线性组合Jacobi迭代与上一次迭代值来实现。
迭代值为Xk+1 = γJacobi(Xk, B) + (1-γ)Xk,其中γ在0-1范围内。
当γ小时,迭代值接近Xk,对应加大权重,收敛慢。当γ大时,迭代值接近Jacobi迭代,收敞快。
放松系数γ决定迭代矩阵R的特征值范围。当γ为0.5时,特征值范围萎缩一半;当γ为0.8时,范围更小。
Multigrid方法将放松Jacobi迭代用于预平滑(pre-smoothing),以平滑近似解。一般采用5步迭代以达到平滑效果。
预平滑后要解错方程在粗网格上,再进行后平滑(post-smoothing)。平滑器的作用是排除近似解中的震荡部分误差。
平滑器处理后,近似解中的震荡误差会被减少。误差方程a∙e=R描述了残留误差。
残留R通过计算B-a∙u得到。这也称残差,可以用来判断收敛性。
将R通过插值下采样到粗网格。粗网格上的R定义了误差方程右侧。
在粗网格上使用多网格迭代法解误差方程,求得e。
将e通过上采样插值回到细网格。近似解u更新为u+e,实现误差校正。
迭代使用相同平滑器后平滑,获取光滑近似解。完成一个多网格循环。
奇数次迭代需判断网格点数是否大于2,满足才下采样到粗网格。下采样和上采样采用平均插值。
多网格算法可以应用于网格点数不是2的指数整数倍时,插值和下采样运算变复杂。
一般多网格将网格细化程度减半,但也可以适度增大网格间隔以减少运算成本。
平滑迭代次数没有一个定值,需视问题和平滑算子而定。不够迭代可能导致误差下采样质量差。
后平滑可以再次应用同一平滑算子,确保最终解平滑。
γ值应使难以在粗网格表达的特征值快速收敛至零。通常取2/3。
测试函数重要,可以发现错误并修正.多网格算法应用一次就能收敛到小误差。
多网格依赖差分方程,无法直接用于线性方程组.但可以使用常用迭代器如互补梯度法.
用户可依问题设置网格细化程度,平滑次数和γ值,以求优化收敛速度。
多网格循环包括细网格预平滑、粗网格求解、细网格解延长和后平滑等步骤。
预平滑降低高频误差组成。求粗网格误差方程,再解延长到细网格增加信息。
多网格递归求解各级误差方程。初次求粗网格误差方程,再递归求更粗网格误差方程。
求解粗网格误差方程后,解延长求出误差,与上层平滑解相加得到次优解,逐层解延长至原问题。
解延长采样误差导致高频部分,后平滑去除高频成分使解更连续。
代码实现采用递归函数,袋各级误差方程,隐式体现多网格循环流程。
测试用例加载右端源项矩阵,从零解起点开始迭代,直至满意程度收敛解决原问题。
讲师设置了多个断点来调试多网格循环中的每一个步骤。当程序运行到断点位置时,会暂停,程序各个变量的当前值会显示在画面左下角。
多网格循环递归地调用粗网解法来求解误差方程。
预平滑:使用平滑算子对当前解进行预平滑,抑制高频错误。
误差计算:计算当前解与真实解之间的误差 Residual。
误差下降:将误差降低到粗网。
粗网求解:在粗网上递归求解误差方程。
解下降:将粗网解下降到细网。
post平滑:对细网加和粗网解后的结果进行后平滑。
循环中每个网格大小都变小,但计算的频率内容越来越低。最后收敛到最细网格上得到好的近似解。
多网格循环每次递归调用时,变量值会变小。随着循环次数增加,解图像中的细节内容越来越低频但清晰度越来越高。
预平滑可以抑制高频错误,错误方程可以提供缺失的高频内容。两者相加获得最好的近似解。
V-网循环只在粗网上调用一次多网格循环,而W-网循环在粗网上调用多网格循环两次。
将程序从V-网循环改为W-网循环只需要增加一次多网格递归调用。这使得每下降一个网格层级,粗网上的多网格循环调用次数增加一倍。
构造1D泊松方程的Jacobi迭代矩阵A。
不同网格尺度下,A的特征值分布不同:
细网格下,特征值接近-1,Jacobi收敛较慢;
粗网格下,特征值接近0,Jacobi能一次收敛。
这是多网格方法在粗网上更快收敛的原因之一。
添加非自然松弛因子γ后,Jacobi迭代的收敛速度为(1-γ)^n,γ越小,收敛越快。
非自然松弛Jacobi迭代可以很快在粗网上收敛低频误差,这也是多网格方法的一个优势。
多网格方法每调用一次,可以减小残差但不会将其减为0。经过多次调用,残差会越来越小。
Jacobi迭代矩阵A的特征值分布与网格大小有关:
细网格下,特征值接近-1,收敛较慢;
粗网格下,特征值接近0,能一次快速收敛。
这是多网格方法在粗网上更快速收敛的原因之一。
加入非自然松弛因子γ后,Jacobi迭代的收敛速度为(1-γ)^n,γ越小,收敛越快。
对泊松方程Jacobi迭代矩阵R的最大特征值与1的差值,随网格大小n的增大呈1/n^2下降。
即使网格越大,收敛速率越慢,网格增大两倍则收敛速率减半。
网格继续粗化到一定程度后,直接求解线性系统比迭代更快。并行计算中也需要在细网转为其他求解方法。
高斯塞德迭代相比雅可比迭代,在解a矩阵时,将对角线和下三角矩阵同时放在左侧,只将上三角矩阵移至右侧。
解一个下三角矩阵其实和解对角矩阵一样简单,即从上到下、从左到右逐行求解。
可以直接将雅可比平滑替换成高斯塞德迭代,前后平滑都使用高斯塞德迭代。
实验结果显示,高斯塞德迭代能比雅可比迭代更快地减小残差。
高斯塞德迭代是非同型迭代,需要按顺序求解每个点,难以并行计算。相比雅可比迭代,计算复杂度可能更低,但并行效率不如雅可比迭代。
平滑器的目的是足够快地减小高频误差分量,如果平滑效果不够好,将会影响多网格方法的整体效率。
当从细网格插值至粗网格时,高频内容可能被错误地表示为低频内容。这称为抽样错误或别名错误。
例如一个细网格中的高频波形,当每隔两个点插值到粗网格后,波形会变成不一致的样子。
如果预平滑未充分减少高频残差,那么在粗网格上求解误差方程会加入不属于原问题的低频项,最终可能增加而不是减少总误差。
当从粗网格回插值到细网格时,线性插值与实际值之间也会有差异,这就是插值错误。插值错误主要属于高频项。
后平滑很好地减小插值错误中的高频项。所以加入后平滑步骤可以帮助多网格方法收敛,尽管从理论上来说后平滑对收敛不是必须的。
∂u/∂t=∂^2u/∂x^2可离散为Du/Dt=a×u ∂^2u/∂t^2=∂^2u/∂x^2可离散为D^2u/Dt^2=a×u
对倾向于stiff问题,采用隐式方法比显式方法更好,如细网格情况下矩阵a条件数大。
后向欧拉法对∂u/∂t=∂^2u/∂x^2离散得到(I/Δt-a)u_{n+1}=u_n/Δt的线性系统。
∂^2u/∂t^2=∂^2u/∂x^2离散得到二阶矩阵方程,经高斯消元后可变为解闹钟矩阵线性系统。各种迭代方法依然适用于该系统。
将隐式方法应用于项目一二中的不定态问题,如后向欧拉法和双步走法。
离散误差定义为解析微分与数值微分近似的差值,即τ=Δ-Δh。
稳定性定义为误差范数与离散误差范数的比率在网格尺寸收缩时保持有界。
隐式方案相比显式方案对stiff问题更好,因为隐式方案具有更好的稳定性。
对不定态PDE,误差在时间上累积。误差方程为en+1=M×en+τn,其中M为离散时间演化矩阵,τ为时间和空间离散误差。
M矩阵特征值范围在-1到1之间时数值方案稳定。特征值大于1时,误差会随时间而越加放大。
对简单推进方程∂u/∂t+∂u/∂x=0进行误差分析,得到误差方程。
分析波动方程∂u/∂t+∂u/∂x=0采用中点差分表示空间微分,前向埃尔利表示时间微分。
假设解在无限域或周期边界下为u(i,n)=Ae^(iωk),ω为频率,k为波数,利用此形式进行分析。
将差分算子施加在正弦波上,解保持正弦波形式但放大倍数α。α称为特征值,正弦波为特征向量。
如果特征值α模值小于1,数值方法是稳定的。采用中点差分和前向埃尔利,α不在单位圆内,方法不稳定。
采用中点差分和后向埃尔利,α在正弦波形式下计算为1+2Δxcosk,其模值大于1,方法稳定。
通过冯诺依曼分析,可以评价不同时空离散方式下的稳定性,如中点-前向不稳定,中点-后向稳定。
对于非线性偏微分方程,稳定性分析需要先将方程线性化。将精确解和数值解的差值定义为误差项。
将精确解和数值解代入非线性算子并相减,得到的是一个带有误差项的线性误差方程。
将误差方程中的非线性算子用泰勒级数展开近似为数值Jacobi算子乘以误差项,要求Jacobi算子的特征值在单位圆内。
线性化处理还会引入一个额外的误差项,但只要Jacobi算子满足稳定条件,这个误差不会积累导致数值解无界。
同样的方法也适用于隐式方法和其他时间离散方式。此外,误差方程中的误差项也可包含由线性化引入的误差。
实际物理模型中,也存在许多对数值方法线性化后的方程本身就是不稳定的模型,例如混沌系统。
混沌系统指的是非线性解本身是合理的,但其线性化后的方程解会发散趋于无穷大。这类系统对常规误差分析不再适用。
对于混沌系统,经验上知道某些不太稳定的数值方案可能更适合模拟其动力特征,但如何分析这类方案的准确性仍是一个开放问题。
模拟物理不稳定模型,特别是混沌系统,需要从新的视角出发重新研究数值方法的分析方法。这将是未来数值PDE研究一个重要方向。
保守定律描述了某些物理量在物理系统中的保守性,如质量、能量和动量等在给定控制体积内的总量不会改变。
标量保守定律中,rho表示单位体积内物理量的密度,是一个标量函数。
保守定律可以用积分方程或偏微分方程表达。掌握高斯定理后,可以将表面积分转变为体积积分,从而转化为偏微分方程形式。
保守定律考虑了流通量F和源项S对系统内物理量的影响。F表示通过控制面流入或流出的物理量,S表示系统内部产生或消耗的物理量。
标量保守定律方程只与ρ有关,ρ是一个标量空间时间函数。求解这个方程可以得到ρ的分布情况。
Finite volume方法是实现数值解保守定律方程的重要方法,可以在结构或非结构网格上针对各种保守定律问题求解。
保守定律可以表示为偏微分方程ρt+∇·f=s。其中ρ表示密度,f表示流量,s表示源项。这称为守恒形式,直接从物理保守定律推导而来。
在一维情况下,∇运算简化为导数运算。 将导数改写为ρ的导数,可以得到原始形式。
在多维情况下,原始形式为ρt+(∇f/∇ρ)·∇ρ=s。其中∇f/∇ρ是矢量。
保守定律也可以表示为ρ的积分形态ρtdx+∫Γf·nds=∫sdx。这是从物理意义出发的初始形式。
ut+ux=0是一维线性标量保守定律方程。其中f(u)=u。
ut+ux=0也是一维标量保守定律,但其中f(u)=u2/2,表示速度为u时的流量。
描述油水在多孔介质中的流动。u表示水体积分数,f(u)表示流量,考虑了复杂动态关系所带来的非线性。
布格斯方程为ut+uux=0。其中u表示解,u也代表速度。方程描述压力波的传播。
解通过有限差分方法数值求解。左边显示实时求解结果,右边显示解在空间时间图中的演化过程。
可见不同区域以不同速度运动,形成不连续性。速度正负区域互相碰撞,形成激波。
激波内分子扩散占主导地位,可将动能转化为内能,因此激波具有消散性,幅值会不断降低。
多个激波形成N型形态,这与超音速飞机产生的声爆现象吻合。
算法长期演化后,解将趋于匀平。不同区域以相同速度运动,不连续性和激波全部衰竭。
方程中的速度相关于解本身,造成各区域运动速度不同。这是布格斯方程中不连续性和激波产生的关键原因。
m代表x方向上的网格点数目或间隔数目。n代表y方向上的网格点数目或间隔数目。
问题1和问题2.1给定明确m和n值。问题2需能接受任意m和n作为输入。m和n一般为2的幂次。
c函数表示解析解,它是x、y、t的函数组合。给不同a、b、c值可以表示不同解析解,如止动波解。
解随时间变化,边界仍按解析解形式变化。非零初值下有波动解,零初值下有止动波解。止动波解需要u0大于1,形成顶部和底部峰值。
给大时间步长,迭代至止动态,解应与给定边界和初值条件下的解析解一致。这可以测试差分场景的正确性。
积分型形式直接描述物质守恒定律。取控制体Ω,控制体内物质量积分对时间求导,等于通过边界δΩ的通量积分加上体内源项积分。
将积分型形式用错断公式转为体积积分。每个空间点,物质密度对时间求导加通量散度等于源项,得到微分方程形式。
将守恒型形式通量项用链式法则变形,得到物质密度对时间求导加通量梯度等于源项。
通量f(u)=u^2/2。满足积分形式、守恒形式和原始形式。
解随时间演化呈直线变形,说明连续区域内导数为零。形成不连续激波处,解呈不连续变化。
对布格斯方程的原始型形式∂u/∂t=u∂u/∂x进行分析。
在(x,t)平面上画一条直线x=ct。沿此线解u是时间t的单变数函数。利用链式法则得到du/dt=c。
对布格斯方程,有c=u。则特征线方程为x=ut+x0,解不变变化。
根据特征线,可以预测解随时间在滑动、斜缩变形。特征线的斜率表示通量函数决定解的运动方向。
利用特征线,可以分析初始值问题的连续解。解保持不变沿特征线传播,形成非连续激波会汇聚。
一般形式为∂u/∂t+∂f/∂u∂u/∂x=0。特征线方程为x=f’(u)t+x0,解不变变化。
连续解可以用特征线分析,当解发生不连续时形成激波。激波区域∂u/∂x不存在,这是解首次不符合微分方程。
原始形式方程在激波处失效,保守形式也失效,因为∂u/∂t在激波移动时不存在。
积分形式是惟一不失效的形式,它将局部不连续区域合并为一个系统整体,时可导数情况下仍保留物质守恒规律。
利用控制体积,在激波两侧解值不连续,通过物质守恒公式可求出激波速度等于左右两侧解平均速度。
例如UL、UR分别表示激波左右解,激波速度us=(UL+UR)/2。初值设定可使激波起初向右后弯向左,表现出弯曲激波运动规律。
学生通过设计不同初值函数,模拟单个弯曲激波的传播轨迹。通过分析初值函数左右侧解值的平均速度,预测激波初期运动方向,并随特征线碰撞而弯曲。
对于一般守恒定律,其通量函数F(u)不再是u^2/2。特征速度改为∂F/∂u。
采用积分形式,在激波两侧近似取微分,利用守恒性得出激波速度公式:
us=(F(uR)-F(uL))/(uR-uL)=∆F/∆u
即通量函数在解值间的差异除以解值间的差异。
在Burger方程特例中,此公式简化为特征速度平均。此结论可推广到任意标量守恒定律。
有限体积法能够准确捕捉解在平滑区域和激波处的行为。而有限差分法仅适用于平滑解,对不连续区域会产生无限大的近似误差。
对于具有可能出现不连续性的非线性守恒定律,选择守恒定律的积分形式。
将求解域划分为多个大小相等的有限体积。体积间界限采用格点。
保存每个体积内解函数的平均值。通过线性扰动近似体积边界解值。将原积分方程转化为微分方程组,利用数值方法求解。
保存体积平均值
重构体积边界解值
计算边界残差
得到离散微分方程组
动差法求解方程组
以上方法可应用于均匀格网和非均匀格网,且在处理激波等不连续性时表现优异。
有限体积法近似求解体积内解函数的平均值,无法保持激波的单值不连续性质,但可以将激波描述为一带状区域。
以边界平均值和左边、右边相邻体积平均值的平均值来近似重构边界解值。
取Δx为体积宽度,求解每个体积内解函数随时间的变化率等于相邻两个体积间数值通量的差。
如果已知边界点解的值,可以直接计算相应边界体积的数值通量,而不用近似重构解值。
考察中心差分方案应用于布格斯方程的离散化过程,得到一个ODYDE系统。但此方案仅适用于平滑解,产生激波时会出现数值振荡。
下次课将讨论如何设计稳定的有限体积方案处理激波问题。
连续区域中解呈直线变化,即特征线;不连续区域中存在激波,即解出现跳变。
特征线速度为F’(u);激波速度为[F(u右)-F(u左)]/(u右-u左)。
原定律可写为F’(u)∂u/∂x=0,称为原型形态,F’(u)即特征线速度。
保存每个计算单元内解函数的体积平均值,近似计数单元界面解值,求解离散形式得到时间方程。
计算单元体积平均值的时间变化率等于相邻单元间数值通量差。
定义计算单元和界面
计算体积平均值
近似求解界面值
求解离散时间方程
以追踪解在不连续区域的传播为目的。
定义计算范围为[0,1],将范围分为32个等距计算单元。计算单元接口点位于网格点的中点位置。
采用正弦函数作为起始解函数分布于计算单元中心点。
设计函数计算时间步长求解的时间变化率du/dt。函数输入参数为解向量u。
本例为Burgers方程,确定单元接口处流量函数表达为u^2/2。
采用固定边界条件,边界点值定为0。
根据中心差分方案计算内部界面流量,边界面流量定为零。
根据有限体积方法,求相关单元界面流量差值,除以单元尺寸,即为时间变化率du/dt。
采用MATLAB内置函数ode45进行时间步长求解,获取数值近似解。
选择时间步长为0.01。
采用for循环,设置5000个时间步长,周期为5秒计算。
构建函数burgers_ddt()返回时间变化率du/dt。
调用MATLAB内置ode45函数,采用burgers_ddt()作为rhs函数,实现时间步长求解。
绘制初始采用正弦函数分布的条件。
根据设置 repeatedly调用ode45进行时间步长求解,绘制迭代结果。
前期计算结果正常,随后形成不连续锐波后,解 funkcion开始明显变形振荡,表明中心差分法无法正确 capturing不连续性。
中心差分法取体积平均值,无法准确描述解在不连续处的行为。
Upwind方法根据风向(扩散方向)选择单侧值,取扩散处的左侧值或右侧值,而不是取体积平均值。
根据f’(u)正负判断扩散方向,f’(u)>0时取左侧值,f’(u)<0时取右侧值。boogers方程中f’(u)=u,直接用u值判断。
采用if-else语句,根据f’(u)值选择取左侧值或右侧值,构建Interfaces_flux矩阵。
采用Upwind方法重新求解,结果表明不连续性保持不变,震波逐渐消散,满足物理特性。
Upwind方法能正确描述不连续性伴随的熵增特性,保留扩散信息,求得的解符合物理意义。
高阶upwind方案适用于波速度一致的情况,如果波速度有较大差异,会产生数值震荡。
多步骤方案可以用一阶显式方法起步,满足时间一致性要求。
Upwind方案仅取单侧的值近似界面通量,会产生顺序误差,导致首主顺序。更高顺序需要引入limiter控制近似。
Limiter可以检测局部最大最小值,在不连续点附近切换为一阶方案,在光顺区域采用高阶方案,实现高低阶混合。
可通过线性化离散方程,研究线性误差方程对误差的进展和衰减规律,也可以求解线性误差方程的伴随方程,分析数值误差对最终量的影响。
高阶Upwind在波速差异大时会产生震荡,limiter可以实现高低阶混合,更好描述不连续性。线性化误差为分析数值误差的关键方法。
布尔格斯方程的解除有不连续点,还可能有多个解,其中扩散波是熵增加的实质解。
特征线表示信息传播方向。从不连续点出发的特征线代表信息产生,符合物理;只有向不连续点汇集的特征线代表信息减少,是非物理解。
Upwind方案近似界面通量时,不连续点正好落在网格接口,界面通量均等,导致时间导数为0,得出非物理 stationary解。
Upwind方案近似错误,产生非物理解。正确应视不连续点为极限解,特征线应originate于不连续点而非terminate于不连续点。
熵条件从物理上决定了特征线应originate于不连续点,但数学上只需满足特征线不能从不连续点产生即可。
布尔格斯方程的物理解没有特征线产生于不连续点。
理想数值近似方案应该能给出物理解。
里曼问题研究不连续点两侧初始条件在无限小时间内的解析解。
下风方面近似界面通量时,可能给出非物理 stationary 解。
古多诺夫数值流量通过解析求解里曼问题得出界面通量,保证得到物理解。
同向传播、反向传播、扩散、激波,求解后确定界面通量。
布尔格斯方程中,激波速度根据左右两侧变量的差决定移动方向。
特征线表示信息传播方向,古多诺夫方案近似特征线走向,得到与物理一致的数值解。
继续研究标量守恒定律,即布尔格斯方程仅一个变量。将扩展到系统守恒定律和二维三维情况。
将解函数存储在每个有限体积内的平均值,而非函数本身。
通过解析求解里曼问题,确定界面通量,获得物理一致解。
同向传播、反向传播、扩散、激波。决定根据左右两边变量和数值梯度的大小关系。
决定由 deltaF/deltaU 的正负号,即流函数 Wyoming 左右两边的差异。
当 characteristics 从负到正,即流函数 Wyoming 曲线最小点,取流函数在此点的值作为界面通量。
将四种情况概括为:(1)特征线同向传播,(2)特征线方向不同,包括扩散和激波。
同向传播:取最大值或最小值
反向传播:取最大值或最小值
激波:根据δF/δU的正负判断,取左值或右值
擴散波:在范围内找到F的最大值点
如果F是单调递增或递减,直接取边界值;否则需考察整个范围内F的极值点
如果δF/δU正,激波向右,取左值;如果δF/δU负,激波向左,取右值
只有FNon单调时才有扩散波,此时取F在范围内的最大值点
如果UL>UR,取F范围内的最大值;如果UL<UR,取F范围内的最小值
如果F非单调,需寻找范围内F的实际极值点作为界面通量
如果F包含导数项,将导数项和非导数项分离,分别处理
根据物理解的左右边界值和特征线方向来确定模拟界面通量数值。
设置时间间隔∆t和空间间隔∆x
初始化解U和时间T
判断每个网格单元内U的左右边界值大小关系
根据汇聚或分散特征线方向,取边界值最大或最小值作为接口通量
用接口通量更新U,计算新的解
重复迭代计算,输出解序列
如果场非单调,需找到实际极值点作为接口通量,而不是简单取边界值
示意扩散波案例,设置Burger方程初值条件。
通过Godunov方案计算数值解序列,并画出解曲线。
可以看出数值解很好地描述了原问题的变化规律。
计算时可能会遇到Python类型错误需要调整。
本次演示结果还需要最后画出解曲线进行验证。
在应用Godunov方案求解波动问题时,设置边界值为0容易导致数值误差。
如果特征线向边界外传播,强加边界值为0会导致数值解错误地取边界值,而不是内部解。
在边界外添加假设单元,使边界索引与内部一致
将边界索引循环范围置为1~n+1,包括假设单元
根据特征线方向,取边界或内部解作为数值通量
修正后数值解在边界处趋向0,并且特征线方向一致,问题得到解决。
Godunov方案是一阶精度,网格细化后精度会提高。但内部可能存在一定跃变。
Godunov方案同样适用于多元难题,只需处理不同方向上的特征线。
数值上可能仍表现为跃变,但由于粘性影响,实际情况可能更加平滑。
在一阶Godunov方案中,我们使用了体素值的平均值近似函数。但是,我们也可以使用更高阶的近似函数。
体素值的平均值实际上是函数中心值的二阶近似。我们可以利用这个特点,在每个体素内用线性函数进行重构。
重构的思想是:使用左右两个体素的值,在接口处进行线性Extrapolation得到重构值。
具体来说,对左侧重构值: 用当前体素值做常数项,用左右体素值差除以网格大小作为系数项,得到线性函数近似。
对右侧重构值同理: 用当前体素值做常数项,用右左体素值差除以网格大小作为系数项,得到线性函数近似。
这样得到的左右重构值,可以直接代入Godunov方案计算界面通量。实现二阶精度。
这个方法利用了体素平均值与中心值关系,以较少的信息实现了二阶重构,可以很好提高求解精度。
在二阶重构方法中,我们需要重构解函数的值ul和ur,分别表示接口左边和右边的重构值。
我们可以使用线性插值法进行重构:
对ul的重构: 使用当前体素值U(I)作为常数项,利用当前体素和前一个体素值差除以网格间距Δx作为系数项,得到线性插值方程。
对ur的重构: 使用当前体素值U(I)作为常数项,利用当前体素和后一个体素值差除以网格间距Δx作为系数项,得到线性插值方程。
重构后得到ul和ur值,就可以直接代入Godunov方案计算界面通量FL和FR了。
为了验证二阶方法,我们实现了Burgers方程和线性波动方程的示意例子。结果表明,二阶方法能很好地处理激波,而没有一阶方法中的不连续现象。
但是,对于线性波动方程,长时间求解会产生无真实性质的振荡。这是需要用限值器等方法来解决的问题。
二维中,守恒定律可表示为:
∂u/∂t + ∂Fx/∂x + ∂Fy/∂y = 0
其中Fx和Fy分别表示沿x和y方向的通量。
将该微分方程应用于有限体法时,需要将其转换为积分形式:
∫∫ω [∂u/∂t] dxdy + ∫∂ω (n·F)ds = 0
这里ω表示控制体积,∂ω表示其边界,n表示边界法向量,F表示通量向量。
通过格网化,控制体积ω可以表示为单元格。单元格内部平均值u称为室平均值。
边界积分可划分为各边上的线积分。我们需要求解每个接口上的平均通量。
我们可以利用临近两单元内部室平均值,建立一维离散问题近似平均通量。
重构法可用于提升精度。但目前我们考虑简单的定值近似。越过接口,看作一维定值问题求解。
这样就将二维问题转化为许多一维子问题,利用一维方法求解各接口通量,从而求解二维问题。
对于二维中有限体法问题,可以有以下讨论:
可以将解进行不同模式的投影获取不同类型的数值方案。常见的有投影到常数模式得到体平均值,或投影到常函数得到节点值。
可以将解投影到多项式模式,如投影到零次多项式得到体平均值,投影到高次多项式可以得到更高精度结果。
对标界面求积分可以采用高阶立面重构或高斯点积分法提高精度。
对标界面求通量可以直接取临近单元体平均值近似,也可以采用重构法利用邻近单元信息进行一维重构得到更高精度结果。
对标界面通量,标量守恒方程可选用阿普温定量法或优等定量法。优等定量法在某些情况下可以避免产生非熵满足解。
可以将单元边 further 分割,或者采用高斯点积分法近似边积分,实现精度的提高。
重构精度高时,可以在界面的多个点进行重构然后采用高斯点积分,此时单元内解需要用更高精度近似。
对未结构网格问题,后续会讨论其中的有限体方法。本次项目考虑结构化网络。
保真定律系统是一类描述物理系统的偏微分方程。它描述物质或者能量在系统中如何保留和流动的定律。
常见的保真定律系统包括:
这些方程都是非线性偏微分方程,包含多个耦合的未知量。比如质量保留定律包含密度ρ作为未知量,而动量保留定律的动量通量则包含ρ和速度u,需要多个方程联合解决。
质量保留定律表示质量在系统内的变化等于流入量减流出量。
用微分形式表示为:
∂ρ/∂t + ∇·(ρu)=0
这里ρ表示密度,u表示流速。
动量保留定律描述系统内动量的变化。系统内动量表示为密度ρ乘以速度u。
动量通量为密度ρ乘以速度u再与自身矩阵积u。加上压力项ρ。
微分形式为:
∂(ρu)/∂t + ∇·(ρu⨂u) + ∇p = 0
这里p表示压力,与ρ、u、e等有联立关系。
能量保留定律描述系统内总能量的变化。总能量表示为密度ρ乘以单位质量内能e。
微分形式为:
∂(ρe)/∂t + ∇·(ρeu + pu) = 0
浅水方程考虑水较浅情况,只有两个未知量水深h和流速u。
微分形式为:
∂h/∂t + ∇·(hu) = 0
∂(hu)/∂t + ∇·(hu⨂u) + ∇(gh2/2) = 0
这些方程都属于非线性偏微分方程系统,难以用解析方法求解,通常需要数值方法求解。
分析一维系统保真定律的代数形式如下:
∂ω/∂t + ∂f(ω)/∂x = 0
其中ω是一个向量,表示系统中多个耦合的变量。f(ω)表示对应每个变量的流动率。
对于标量情况,可以通过求解微分方程∂f/∂ω=0得到标量特征线方程。
但对于向量方程,求解∂f/∂ω时会得到一个雅可比矩阵。
如果雅可比矩阵为对角矩阵,表示每个变量独立,方程可分解为多个标量方程。此时特征线速度等于对角线元素值。
对一般非对角矩阵,需要进行特征值分解。对原方程乘以特征向量矩阵u,通过特征值Λ可将方程对角化,得出一个行波速度不同的耦合标量方程组。
特征线给出了解的传播方向和速度。但与标量情况不同,多变量系统特征线可能会互相影响,生成复杂的转换。
只能通过数值方法求解多变量系统,特征线信息对构建数值差分方案至关重要。
为多个耦合变量构造界面流量是一个复杂问题。除特殊情况外,通常需要求解瑞利问题来获得交界面流量。
考虑添加一个数值消散项到系统中,使其近似物理弱解。数值消散起到类似物理消散的作用,但不需消散趋于0极限。
将界面流量近似为两个点流量平均值加上数值消散项。数值消散量B取最大速度谱值或特征值绝对值。这样可以获得上风差分性质。
数值消散项相当于加入一个差分项,发散量与网格间隔成正比。而最大发散量B仅与方程有关。随网格间隔缩小,发散量趋于0,近似了物理弱解。
| 如果流量F与变量W成比例,系数为最大风速u,则B= | u | ;一般情况下,B取F对W求导数Jacobi矩阵特征值绝对值最大值。 |
数值消散法解决瑞利问题,给出界面流量和上风性质。随网格细化,趋于正确物理解。
W表示保存量,单位为X;F表示通过单位面积单位时间保存量流量,单位为X×速度。
如果F与W非线性关系,应考虑∂F/∂W的单位为X/时间,即速度。B的单位也是速度。
如果有真实粘滞项,总流量公式为:
F(W) - μ × ∂W/∂X
其中,∂W/∂X使用二阶有限差分近似。
将粘滞项∂F/∂X离散后,保持方程的保质量。但粘滞项内部的一阶导数使用有限差分近似。
这样即包含了数值粘滞和物理粘滞,同时保持方程数值保质量。
在前面讨论布格尔方程例子时,作者忽略了设置周期边界条件时应复制数组端点的值。这导致数字场中出现不物理的振荡。
1单变量保存方程的总体变化定义为解在空间整体变化的积分值。
2单变量保存方程解随时间的总体变化应不增长,并在有震波时单调减少。
3如果数值方案不能保持数值解在时间上总体变化不增长的性质,会产生不物理的数值振荡。
如果线性单变量保存方程的数值方案是线性的,那么要求总体变化减少(TVD),它的顺序只能是第一顺序。
但如果方案是非线性的,即用非线性差分方程离散线性微分方程,则可以高阶而保持TVD性质。这需要利用限制器等非线性机制。
考虑Partical U/Partial T + Partical U/Partial X = 0这一简单线性偏微分方程。用有限体积法离散后,可以得到无限条特征线。假设初始条件有不连续点,则随时间特征线沿着不变形地前进。
在有限体积法中,进行重构后可能会超出初始解的最大最小值,产生不物理振荡。例如在局部极值点附近进行线性插值重构,可能会导致重构值超出实际范围。
限流器是一种非线性机制,用来检测局部极值点。如果发现局部极值,则采用常数重构;如果发现区域内解为线性分布,则采用线性插值重构。通过限流器,可以在光滑区域使用高度准确的插值,在不光滑区域使用稳定的重构。
限流器的作用是让重构方式在不同区域连续地进行切换。避免重构方式的突变会损害数值解的准确性。实现限流器需要考虑如何检测局部极值、如何在不同重构方式间进行平滑切换等问题。
限流器是根据UI和相邻单元的值计算出重构梯度的一种非线性方法。它可以检测局部最大值和最小值,在光滑区域使用线性插值,在不光滑区域使用常数重构。
重构梯度根据UI,UI-1,UI+1的值计算:
如果UI在UI-1和UI+1之间,表示在光滑区域,则使用(UI+1-UI)/(UI-UI-1)计算梯度;
如果UI为局部最大值,则梯度为0;
计算梯度同时考虑不超过相邻单元的值。
限流器函数φ取决于比例R=(UI+1-UI)/(UI-UI-1) ,它必须满足:
这些约束定义了φ在R坐标轴上允许的区域。
使用计算出的梯度和限流器函数φ,可以根据UI-1,UI,UI+1的值重构出UI周围单元的值:
梯度 = (UI - UI-1)/(Δx) * φ(R)
从而实现在不同区域动态切换线性插值和常数重构,防止产生非物理振荡。
限流器函数φ(R)定义为:
它采取最保守的重构方式,实现最小斜率重构。
限流器函数φ(R)定义为:
它采取最激进的重构方式,实现最大斜率重构。
限流器函数φ(R)定义为:
它在最小模式和超限流器之间找一个平衡,在0-1之间以光滑曲线连接。
所有限流器函数都满足φ(1)=1的约束,以实现二阶数值差分约束。它们会在光滑与不光滑区域自动切换重构方式。
有结构网格在处理复杂几何形状时难度大,但采用无结构网格可以很好适应复杂几何形状。
网格结点的形状和位置不受限制,可以适应任意几何形状。
无需计算微分算子,只需要计算单元界面处的流量即可。
可以使用网格细化,在解算精度高的区域增加网格密度。
网格交错重叠区域可以描述复杂体积形状。
三角形和四边形单元可以适应 near 边界区域的曲面。
网格密度可以根据解算精度的需要进行自适应改变。
单元平均值等于单元所有界面处流量的和。对于具有多个界面的单元,其贡献等于所有界面流量之和。
这种方法可以在任意无结构网格上求解守恒型微分方程。
PD工具箱可以使用图形界面绘制网格,如圆形、矩形等,也可以设置子网格名称和边界。可以导出网格信息,如点集P、边集E、三角形集T等。
Draw点集P、边集E、三角形集T后,可使用工具提取内部边,从而得到全部边集。
每条边的法向量由端点位置计算得出。法向量方向为从第一个端点指向第二个端点的反方向。
绘制点集P,绘制边集E以连接点,绘制法向量集以箭头表示法向。
1 计算每个单元的单元切片率
2 循环计算每个界面的切片率
3 将同一切片率同时添加至一个单元内,从另一个单元内减去
以保证数值连servativeness
4 对所有界面重复复1~3步,算法总体保持质量守恒
5 边界单元添加或减去切片率时,只作用于内侧
uL接口值计算为边界单元左侧界面值,uR接口值计算为边界单元右侧界面值。
通过邻近单元值计算r系数,以适应光滑结构与非光滑结构。
根据r系数值范围应用限流器函数,限制r在0-1范围内,避免产生非物理波动。
如果r大于0,r=r;如果r小于1,r=r;否则r=1。
将限流器函数封装为独立函数linear(r),方便重复调用。计算uL和uR时,都调用该函数,按对称性计算不同方向的r系数。
讲师表示会在下周课上测试修改后的代码,验证是否实现了二阶实质精度的有限体积法。
限流器是引入非线性的机制,即便求解线性方程。通过重构过程限制解在特定区域的取值。
重构把体积平均值推广到接口,但只在解连续光滑的区域适用。在极值点需要抑制重构。
R比为两个邻近差值的比率,反映解的光滑程度。限流器函数根据R比控制重构系数φ的值范围。
如果R≥0,φ=R;如果R<1,φ=R;否则φ=1。
学生指出计算R比的公式错误,应为(u(i+1)-u(i))/(u(i)-u(i-1))。
只需要修改限流器函数即可实现不同限流器,例如VanLeer限流器为φ=(2R)/(R+1)。
若R为无限或NaN,限流器函数应返回0,以避免φ出现异常值。
运用不同限流器求解,验证其TVD特性,如VanLeer限流器产生更尖锐的数值震荡。
重构接口值时,重构左边界和重构右边界分别对应不同的方法。
两种方法只关注解在特定点或体积中的值,而不是构建整个解函数。
若边界条件给定,直接采用边界值。
若特征向外,使用内部点值计算边界流量,无需重构。
若特征向内,可以利用边界值和第一内点构成斜率进行重构。
计算相邻两点之间解的变化率,反映解在该点的光滑程度。
重构时必须区分特征方向,避免强加不合理的边界条件值对内部解造成影响。
回答了有限差分和有限体积方法的异同点,以及重构边界和内部点的不同方法。
有限元素从偏微分方程的弱形式入手进行离散化,而非不同形式或积分形式。
投影是有限元素的核心概念。有限元素实际上就是对函数的一种投影 approximation。
有限差分记忆点值,有限体积记忆体积平均值,而有限元素采用函数投影的方法对函数进行逼近。
连续 piecewise 线性方法,逼近函数求最小积分二乘距离加权和。
非连续 piecewise 线性方法,同样逼近函数但允许在边界处出现不连续。
通过构建套弹性弦的最小能量逼近问题,得到piecewise线性函数空间下对原函数的最佳逼近,从而实现函数的投影。
介绍了有限元素的基本概念和区别,并通过交互示范解释其基于最小二乘投影的数值逼近原理。
有限元素不是直接在差分形式或积分形式下离散化偏微分方程,而是基于其弱形式进行离散化。
投影是有限元素的关键概念。有限元素实际上是对函数的一种投影approxiamtion。
有限差分记忆点值,有限元素记忆函数在格点的取值。有限元素采用基函数进行函数空间的投影逼近。
有限元素求解函数FH(x),使其在基函数空间XH内使函数F(x)与FH(x)积分差的二次型最小。这实现了函数F(x)在XH上的最佳逼近,也即投影。
基函数空间XH具有可能的基,使其任意函数可以表示为基函数的线性组合,从而使XH变为有限维空间,可以在计算机中表示。
举例说明了区分单元内多于一个点的基函数构造方法,证明任意连续平滑函数都可以表示为这些基函数的线性组合。
有限元素的基函数空间可以有不同选择,如平面片断状可不连续的情况下,基函数需要重新构造。
举例说明如果允许函数不连续,可以用三个区间且每个区间内线性,但在区间边界处值可不连续的基函数。这种情况下基函数空间维数为2倍于区间数。
基函数空间不唯一,也可以采用多项式基函数空间,或 Fourier基函数空间等。不同基函数决定了空间维数的不同。
在二三维情况下,可以采用三角形或四面体有限元,把问题分解在每个单元上,并采用单元内线性基函数得到解。这比结构网格更直观。
有限元素方法在选择基函数空间上很灵活,可以满足不同问题需要。基函数选择决定了求解过程和精度,需要根据实际问题特点进行合理选择。
有限元素投影方程可以表示为求解一个二次最小化问题,即使函数在有限维空间内近似值与原函数误差最小。
二次最小化问题的解需要满足由任意 virtual disturbance Delta G引入的函数误差大于或等于0。利用Delta G可以很小的特点,可以推导出投影方程的等式形式。
投影方程的等式形式为:抗扰动Delta G与近似误差之间的内积积分为0。即投影后的近似解与原函数之间的误差向量与空间中任意向量都夹角为90度。
内积为0等价于投影误差向量与任意向量之间垂直。也就是说,投影后的近似解到原函数的连线在有限维空间内是不前进的。这与实际投影效果一致。
利用二次最小化问题的性质,提供了投影方程的等式形式表示方法,它从数学上对应投影误差的正交性,更好地解释了投影的物理意义。
投影方程要求积分残差与任意扰动函数的内积为0。这意味着残差向量与任意向量之间垂直。
可通过积分型换得公式将二阶导数项转换成边界项与一阶导数项,从而消去二阶导数。
投影方程实际上是一个无限个方程组,但如果采用有限个基函数试求解,则只需考虑这些基函数对应的方程。
将未知函数表示为基函数线性组合,将投影方程分别对基函数求积分,得到矩阵方程。
以上过程giving了泊松方程在有限元素方法中的解法。
投影方程期望残差与各基函数的内积等于零,但不一定能完全满足,这取决于选择的近似空间。
如果没有指定边界条件,直接解决偏差方程,矩阵A很可能为奇异矩阵,这是因为问题本身没有唯一解。
通过指定边界条件,例如Dirichlet边界条件u(0)=0,可限制解空间的维数,使问题具有唯一解,此时矩阵A不再奇异。
近似解与投影应使用相同基函数空间,否则方程数量与未知数可能不匹配,矩阵可能为奇异。如果选择不同基函数,必须确保两种基函数生成空间维数相同。
下节课将介绍有限元素方法的更广泛概念,不限于泊松方程,给出其在其他分离型方程中的应用。
有限元方法将无限维函数空间投影到有限维 trial 函数空间中,近似求解泊松方程。
Trial 函数空间采用三维线性元素划分成网格,在每个三角形内部函数线性。函数在边界取值为 0 满足边界条件。
将残差函数投影到 Trial 函数空间中,要求残差内积与任何 Trial 函数都为 0。这相当于将残差函数在有限维空间中的投影为 0。
使用格林公式将求导移到 Trial 函数上,使 Trial 函数只需要一次可导。然后利用积分进行残差与 Trial 函数的内积求解,得到线性方程组。
给出一个三角形单元,定义不同控制节点函数值的线性函数样式,并给出泊松方程在此单元中的有限元形式。
目前采用的只是一次线性 Trial 函数。如果使用高阶元素,需要确保函数和其导数在单元与单元之间能够连续。
利用格林公式将泊松方程中的求导转移到试函数上,使求导形式放松成为弱形式。
将计算区域划分成三角形单元,在每个单元内定义基函数为线性函数。基函数的值在节点为1,其他区域为0,组合成基函数集。
将解函数U和试函数V分别用基函数展开为线性组合。将弱形式中U和V分别替换,得到线性组合中每个基函数对应的 弱形式积分等式。
将每个基函数对应的弱形式积分等式提取出来,构成矩阵方程。矩阵左端代表积分项,向量代表源项积分。求解线性方程组可以得到解函数的配比系数。
根据基函数展开式和各配比系数,求得解函数在每个三角单元内的具体表达式,完成有限元方法的计算过程。
利用随机数生成1D空间内网格点坐标集合X,网格点个数设为8个。将X从小到大进行排序。
基函数定义为在每个单元内以线性函数变化,在节点值为1,其他区域为0。利用X和单位矩阵可以表示每个基函数。
K矩阵用于保存每对基函数梯度内积。分析基函数不相交区域内积为0,相邻单元内积表达式为-1/H。
在1D情况下,源项取常数1,B向量元素为对应的基函数整域积分值,可表示为H1与H2之和。
利用矩阵K和向量B,构建线性方程Ku=B,利用 Matlab命令K\B求解,得到解函数u在每个节点的值。
将计算结果u与X作图,并绘制解析解函数-x^2(1-x)与之对比。结果表明数值解能较好地近似解析解。
在这个视频中,教授们介绍了如何使用Matlab构建 Poisson 方程的有限元素法模拟,完成2D情况下矩阵K和载入向量B的构建。
首先使用Delaunay三角剖分生成一个238个三角形元素的网格,其中p存储节点坐标,t存储三角形节点索引。
然后定义基函数,在每个三角形内部线性变化,其余为0。计算每个三角形内每个基函数的梯度,通过求逆运算得到基函数在三个节点上的梯度值。
构建 zeros(238,238) 的K矩阵储存每个节点对的积分。遍历每个三角形,计算相应三个节点的梯度,通过积分得出K矩阵中的块增量值。
但是K矩阵大小应为节点数而不是元素数。此外,需要考虑边界条件,初始化时包含边界节点,后将其删除。
三角形内积分计算需要考虑三角形面积。构建过程中需要检查矩阵乘法顺序是否正确。
这个视频主要关注有限元法中矩阵K和载入向量B的构建算法以及实现难点,为后续Poisson方程数值模拟奠定基础。
教师首先计算每个三角形元素的面积。利用三个节点坐标组成的矩阵求解行列式,即可得到面积值。
然后初始化大小为节点个数的K矩阵和载入向量B。K矩阵反映节点对之间的交互作用,B向量记录每个节点所承受的载荷。
要逐个三角形元素填充K矩阵。对每个元素,取其三个节点,计算对应的基函数梯度。利用梯度内积和元素面积,计算该元素对K矩阵中对应节点对的贡献值。
由于对易性,只需计算主对角线以上部分。对每个节点对,通过累加各元素的贡献值,逐步填充完整K矩阵。
B向量的计算方式类似,也是累加每个元素面积的1/3份。乘以单一载荷函数得值,即为该元素对B向量的贡献。
循环结束后,即完成了K矩阵和B向量的构建。K矩阵反映不同节点间的耦合关系,B向量记录每个节点所受的整体载荷,为后续求解提供数值模拟必需的数据。
函数空间是向量空间,必须满足向量空间的四大 axiom:
加法交换律
可视为 scaler 乘法
标量与向量的 distributive property
存在零向量
教师给出两个常见函数空间:
L2空间:函数和其差毕函数的平方可无限积分收敛。
H1空间:函数和其一阶导函数的平方可无限积分收敛。
举例说明:1/x 不属于 L2空间,但属于 H1空间。sqrt(x) belonging L2但不属于H1。
高阶函数空间Hk需要函数连续可微至k阶,且k阶导函数平方积分收敛。
函数空间下定义范数来衡量函数大小:
L2范数为函数平方积分开根号。
H1范数为函数和一阶导函数平方和开根号。
给出直观解释:L2下紧凑函数范数小,H1下增大导数会导致范数增大。
常见空间关系:Hk包含于Hk-1包含于L2。δ函数不属于任一空间。
函数是将一个元素映射到另一个元素的关系。对空间中的函数,如果满足线性性质,即:
如果对函数f进行α倍缩放,对应的函数值也会缩放α倍;
如果函数f是g和h的和,那么函数值也等于g和h函数值的和。
那么这个函数就称为线性函数。
函数态是定义在函数空间上的函数,它将一个函数映射到实数。满足以下性质的函数态称为线性函数态:
对函数f进行α倍缩放,函数态值也会缩放α倍;
如果函数f是g和h的和,那么函数态值也等于g和h函数态值的和。
例如对函数求积分就是一个线性函数态。
二线性函数态将两个函数映射到一个实数,它对第一个函数和第二个函数都满足线性性质。
具体来说,如果函数态记为B(f,g),则:
就是说对第一个和第二个变量都满足线性性质。
在某个函数空间中,我们可以定义函数的规范来衡量函数的“大小”。对应地,对于线性运算符L,其规范定义为:
‖L‖=sup‖Lf‖/‖f‖
这里f取值为该函数空间中的所有非零函数,‖f‖为f在该空间中的规范,‖Lf‖为L作用在f上的规范。
定义一个函数空间X,例如H1(0,1)表示0-1区间内的H1空间。
给出一个有界线性形式l和有界双线性形式b,例如:
l(f) = ∫01 f(x)dx
b(f,g) = ∫01 df/dx dg/dx
用线性形式l和双线性形式b定义PDE的弱形式:
b(f,g) = l(f)
对任意f∈X都成立。
对b(f,g)进行积分变换,可以得到原PDE:
d2g/dx2 + 1 = 0
积分变换会产生边界项,要求边界项为0,就得到自然边界条件:
| dg/dx | x=0 = dg/dx | x=1 = 0 |
限定函数空间X不包含某些函数,也可以强制边界条件成立,这就是本本质边界条件。
总之,弱形式定义解除了原PDE中的微分概念,但包含所有物
常规的自然边界条件包括Newman边界条件和Robin边界条件。
Newman边界条件意味着解域的法向导数等于某个定值。
在弱形式方程中,双线形形式不需要修改,仍为整域内积。但是线性形式需要在边界上增加一项。
具体来说,线性形式中的体积积分项不变,但需要在边界上增加一项等于法向导数乘以测试函数v的积分。
这样修改后,弱形式方程在体内满足偏微分方程,同时在边界上也满足Newman边界条件。
Robin边界条件意味着解域法向导数加上一个线性函数等于某一定值。
弱形式方程中,双线形形式不需要修改。线性形式体积积分项也不变。
但是需要在边界上增加两项:一个为测试函数v乘以一个定值的边界积分,另一个为测试函数v乘以解域u的边界积分。
这样修改后弱形式方程在体内和边界上都能满足Robin边界条件形式。
重要的是,两边弱形式方程中的测试函数v都可以取任意函数,这就保证了弱形式方程能满足原始边界值问题。
如果边界条件为非零形式,如u=ub在边界,这时无法通过修改弱形式方程的双线性形式或线性形式来满足这样的边界条件。
我们需要将解空间projec到一个仿射空间中,这样解空间不再是一个线性空间。
具体来说:
定义一个更高维的函数空间XH,解u属于这个空间
定义一个附带非零边界条件的仿射空间XHB,该空间是一个子空间
将XH中的函数f投影到XHB中,投影后的函数记为fH
fH是距离f在XHB中的最短元素
f-fH与XHB中任意元素u-fH正交
这就定义了投影,也构成了弱形式问题中的测试函数空间
弱形式方程设置为,双线性形式不变,线性形式内积项不变,测试函数取自XH0空间
这样修改后,弱形式方程在仿射空间XHB中满足本质边界条件。
弱形式包含双线性形式B(u,v)和线性形式L(v)。
为了进行有限元分析,需要选择两个子空间:
解空间XH,解u属于这个有限维空间
测试函数空间XHB,测试函数v属于这个有限维空间
保持B(u,v)和L(v)不变,将它们应用于有限维空间XH和XHB。
这样弱形式仅需要对XH基进行求解,得到一个有限维线性方程组。
解u在XH中的表示为基函数的线性组合,所以问题转换为求解该线性组合的系数。
选择不同的XH和XHB,可以得到不同种类的有限元方法。例如选择一次连续函数就得到一次有限元,选择不连续函数就得到DG方法。
在实现的时候,需要为各种边值问题详细阐述具体算法,包括稳态问题和未稳定问题。也需要讨论选择不同的函数空间。
弱形式中的线性形式L(v)需要对各种情况下的f(x)进行积分。
将整体积分分解为各元素内部的积分求和,这样f(x)和v(x)在每个元素内都是光滑的。
利用高斯积分法可以很好地近似各元素内部光滑函数的积分。
高斯积分法用预设的立即点和权重近似积分,对多项式以及光滑函数具有很高准确度。
它可以专门逼近n阶多项式,并且随着点数的增加,准确度不断提高。
通过matlab函数可以获取不同数量立即点和权重。
应用高斯积分法计算正弦函数在区间[0,pi]内的积分,结果精度随点数增加而提高。
将弱形式中的线性项积分转换为立即点求和,实现了弱形式离散。这是有限元方法中常用的数值积分技术。
对于非定态微分方程来说,使用有限元方法时,时间导数项将不再是单纯的。施加了有限元素后,时间导数项中将出现一个称为“质量矩阵”的项。
将微分方程投影到一个函数空间中,使左端项减去右端项在该空间中的内积为0。
将空间限制到有限元素子空间中,使方程只在这个子空间内成立。
将解的近似表达通过基函数来表示,将方程离散为代数方程组。
对时间导数项使用离散方法,如欧拉法或龙格-库塔法进行时间积分。
质量矩阵Mij定义为基函数φi和φj在空间上的内积。它是一个对称矩阵。
由于存在质量矩阵,采用隐式时间积分方法较为合理,避免每次时间步幅都需要对质量矩阵求逆。
比如使用梯形法则:
M(Un+1 - Un)/Δt = (1/2)(AUn+1 + AUn) + Bn
移动项集中到Un+1左边,形成一个线性方程需要解。
对于非定态问题,有限元素方法将产生质量矩阵项。这给时间离散带来一定困难。但采用隐式时间积分方法,通过线性方程组的形式解决问题,效率较高。
高阶元素指在单个网格中的函数空间,使用高于一阶(线性)的多项式来逼近求解函数。
例如:
使用高阶元素比线性元素更能真实地表达平滑函数。
二阶元素相比线性元素,可以提高近似精度一个阶,将误差项从$h^2$降低到$h^3$。
这是因为二次多项式能更好地逼近曲线,而线性函数只能实现最佳逼近时为通过中点的直线。
对于单个元素,二次多项式需要3个基函数才能构成完整基:
对于多个元素,需要考虑元素间的连续性。
可以通过线性组合原有的线性基函数,得到二次连续基函数:
这些新基函数与原基函数的线性组合,就可以构成二次连续多元素函数空间的基。
使用高阶元素后,弱形式和实现过程保持不变:
只是基函数从一阶升级到二阶,但原理完全相同。
使用二阶多项式的元素,每个元素需要3个基函数:
对每个元素,求解微分方程需要积分各基函数的导函数乘积。
这个积分可以分解为每个元素的积分求和。
由此可构造全局系数矩阵A,其元素来源于每个单元的local矩阵A^i的叠加。
每个元素的local矩阵A^i都是3×3的。
其元素根据对应基函数导函数的积分计算得出:
如中间基函数导函数的积分需要画图求解。
采用符号计算或高阶定积分可以有效求出各项积分表达式。
高斯定积法可以用少量点高精度求解各基函数在元素内的积分。
这可以有效构造并求解高阶元素问题,且适用于2D和3D情况。
如果微分方程中的系数根据解而变化,就形成非线性偏微分方程。
非线性方程不满足线性叠加原理,难以求解。
利用积分与部分积分将微分方程转换为弱形式方程,以适应较低的解空间光滑程度。
将解空间和测试函数空间限制在有限维子空间内,从而将弱形式方程转换为有限个等式。
用有限个基函数展开解空间中的解,从而将问题转化为求未知系数。
非线性项的系数根据解而变,导致得到的方程组为非线性方程组。
常用的迭代求解方法有牛顿法和平方残差法等,通过迭代逼近得到解。
将非线性方程组写为F(a)=0形式,用雅可比矩阵J=∂F/∂a描述方程组。
利用迭代求解:ak+1=ak-J-1F(ak),初始值a0要足够接近真解。
若求解非线性动态方程,前一时刻解可以作为下一个时刻的初值。
近似计算雅可比矩阵J,ak+1=ak-JapproxF(ak),Japprox需要计算简单。
收敛速度可能比牛顿法慢,但计算复杂度较低。
适用于雅可比矩阵难计算但近似值可取的情况。
考虑边界值为0的单变量无限泊松方程:
∂u/∂x = f(x)
其中f(x) = sin(πx), x属于[0,1]
有限元解会使解的能量误差取得最小值,即梯度误差在整个域上的积分最小。
对于泊松方程,能量误差定义为解梯度在整个域上的面积误差。有限元解能很好地匹配这个误差定义。
很多偏微分方程描述系统处于平衡态,这种平衡态通常对应于系统能量或动量的极小值。
泊松方程描述压强场,其能量函数为解导数平方的积分。解析解使此能量函数达到极小值。
梁方程描述构件变形,其能量与构件弯曲有关,对应二阶导数平方项。
有限元解会使解对应能量函数的能量误差取得最小值。
将偏微分方程写成对称正定弱形式:
a(u,v) + L(v) = 0
其中a是对称双线性形式,L是线性形式。
当a对称正定时,对应能量范数存在。
有限元解会使a范数取得最小值,从而匹配误差定义。
对于泊松方程,能量范数定义为解在整个域上的梯度误差面积。
有限元素方法寻求最小化该能量范数,使解匹配误差定义。
伴随偏微分方程的弱形式满足对称正定性。
精确解和有限元素解都满足该弱形式定义。
对精确解和测试函数,双线性形式的值是相同的。
根据双线性性,有限元素解使误差函数与任意测试函数的双线性形式值等于零。
再根据正定性,任何其他解的能量范数都不可能小于有限元素解。
通过假设另外一个解,展开其能量范数,可得出结论。
其中利用了双线性形式与正定性的性质。
有限元素解使每个元素内解的平均梯度等于元素内精确解梯度的平均值。
这满足最小能量原理,使解在每个元素内都匹配精确解的梯度特征。
解析解和有限元素解相比,通过改变形状可能得到更小的L2范数。
对任意函数f,若在某点或者平均值为0,则其L2范数受其导数L2范数约束。
可以用常数a使能量范数小于或者等于a平方加1乘以H1范数。
通过能量范数和H1范数的相当性,可以证明有限元素解的H1范数误差受常数控制。
L2范数和能量范数不具有相当性,但L2范数误差受H1范数误差控制。所以L2范数误差也受常数控制,但常数与最小L2范数误差的比例没有限定。
有限元素解在H1范数和L2范数下都受常数控制,但仅在H1范数下表现为近似优化。
我们证明了有限元素解在有限维线性空间中实现最小能量范数误差。
可以证明有限元素解在索伯列夫范数和L2范数下都有常数约束。
如果功能G(U)是解U的线性函数,那么我们关心G(U)的误差而不是解U本身的误差。
为估计G(U)的误差,我们定义另一个离散问题——伴随方程。它使用相同的双线性形式,但右端为G(ψ)。
解决伴随方程可以得到伴随解φ。利用φ,可以将G(Uh)-G(U)表示为原问题的残差。
利用有限元素解Uh代入原问题,与φ作为试函数,残差就等于G(Uh)-G(U)的估计值。
如果使用同样的有限维空间解伴随方程,残差为零,无法评估误差。所以需要使用细致一级的空间。
伴随方程不仅可以对线性函数估计误差,后面会讨论它的其他应用。
通过解偏微分方程描述空气动力学,可模拟超音速飞机的阻力和气流噪声。根据微分方程解的灵敏度衍生,优化飞机的几何形状以降低阻力和噪音。
通过解偏微分方程来模拟和分析航天器产生的低频噪音,指导修改设计以减少结构损伤。
理论上可设计高速飞机发动机,但细节需要通过模拟来优化,实验成本高。
通过模拟油气场中的流体动力学,研究不同钻井位置和采油方法,选择能最大限度提高采出率的方案。
还包括模拟飞机、航天器等外体表面流形结构,化工过程等,通过解计算机模拟的微分方程,优化设计和控制变量以满足性能与安全要求。所有这些应用都受益于求解微分方程和计算灵敏度导数。
利用有限差分或有限元等方法,可以将偏微分方程离散化为线性系统Ax=b。
设计优化中的优化目标j通常是线性函数c^T×,优化任务是计算j对设计参数s求导数∂j/∂s。
直接求解∂x/∂s,然后计算c^T×∂x/∂s,适用于s低维j高维场景。
求解伴随方程A^Tφ=c,则有j=b^Tφ,∂j/∂s=φ^T×∂b/∂s,无需求解x。适用于s高维j低维场景。
求解伴随方程只需求解一次线性系统,即使s类别多,也可以一次性计算出所有∂j/∂s。
直接法适合s低维,灵敏式法适合s高维。两种方法都能高效计算出线性系统中的设计参数对目标函数的导数,用于参数优化。
利用函数f(x,s)=0来表示非线性系统,其中x代表变量,s代表设计参数。
如果对x求导的雅可比矩阵非奇异,则在一定范围内存在唯一解x与s的映射关系。
目标函数j可能是x的非线性函数。
通过扩增目标函数j加上λ*f(x,s),得到新目标函数J(x,λ,s)。
∂f/∂xδx + ∂f/∂sδs = 0
∂j/∂x + λ*(∂f/∂x)=0
∂j/∂s = λ*∂f/∂s
只需求解一次非线性方程及一个拉格朗日乘数方程,即可高效计算任意设计参数下目标函数导数,用于参数优化。
利用∇·(a∇u)-f=0表示非线性系统,其中u为变量,a为设计参数矩阵。
目标函数j为场u的非线性函数egral(c·u)dx。
通过扩增目标函数加上λ*(∇·(a∇u)-f)得到新目标函数J(u,λ,a)。
∇·(δa∇u)+δa·∇u=0
∇·(a∇φ)+c=0
∂j/∂a=∫(∇φ)·(∇u)dx,φ为拉格朗日乘数场
通过求解系统及拉格朗日乘数方程一次,便可高效计算任意a下目标函数导数,用于参数优化。