第五讲——显式差分和隐式差分(5)(格式整齐)
显式有限元和隐式有限元
按照计算每一时刻动力反应是否需要求解线性方程组，可将直接积分法分为隐式积分方法和显式积分方法两类。
隐式积分法是根据当前时刻及前几时刻体系的动力反应值建立以下一时刻动力反应值为未知量的线性方程组，通过求解方程组确定下一时刻动力反应。
隐式方法的研究和应用由来已久，常用的方法有线性加速度法、常平均加速度法、New mark方法、Wilson-θ法、Houbolt 方法等。
显式积分法可由当前时刻及前几时刻的体系动力反应值直接外推下一时刻的动力反应值，不需要求解线性方程组，实现了时间离散的解耦。
解方程组一般占整个有限元求解程序耗时的70％左右，因此，这一解耦技术对计算量的节省是可观的。
隐式方法大部分是无条件稳定的，显式方法为条件稳定。
显式方法的稳定性可以按满足精度要求的空间步距确定满足数值积分稳定性要求的时问步距来实现。
显式方法受条件稳定的限制，时间积分步长将取得较小，但计算经验表明，对于一些自由度数巨大且介质呈非线性的问题，显式法比隐式法所需的计算量要小得多。
因此，随着所考虑问题复杂性的增加，显式积分法得到重视。
对于显式与隐式有限元的理解关键字: 有限元显式隐式显式算法和隐式算法，有时也称为显式解法和隐式解法，是计算力学中常见的两个概念，但是它们并没有普遍认可的定义，下面只是我的一些个人理解。
一、两种算法的比较1、显式算法基于动力学方程，因此无需迭代；而静态隐式算法基于虚功原理，一般需要迭代计算。
显式算法，最大优点是有较好的稳定性。
动态显式算法采用动力学方程的一些差分格式（如广泛使用的中心差分法、线性加速度法、Newmark法和wilson法等），不用直接求解切线刚度，不需要进行平衡迭代，计算速度快，时间步长只要取的足够小，一般不存在收敛性问题。
因此需要的内存也比隐式算法要少。
并且数值计算过程可以很容易地进行并行计算，程序编制也相对简单。
但显式算法要求质量矩阵为对角矩阵，而且只有在单元积分点计算尽可能少时速度优势才能发挥, 因而往往采用减缩积分方法，容易激发沙漏模式，影响应力和应变的计算精度。
古典显式差分格式matlab
古典显式差分格式matlab古典显式差分是一种常用的数值方法，用于求解偏微分方程。
在Matlab中，可以使用该方法来解决一维或二维的偏微分方程问题。
本文将介绍古典显式差分格式的原理、步骤以及在Matlab中的具体实现。
1.古典显式差分格式原理古典显式差分格式是一种基于近似的方法，通过将偏微分方程中的导数用差分形式取代，从而将偏微分方程转化为一个离散的代数方程组。
该方法的基本原理是将空间域和时间域划分为离散网格，在网格上近似求解原方程。
2.步骤古典显式差分格式的求解步骤如下：(1)网格划分：将求解区域划分为N个小网格，形成离散的空间域。
(2)时间步长确定：确定离散的时间步长Δt。
(3)差分方程建立：根据偏微分方程，将导数项用差分近似表示。
(4)边界条件处理：确定边界条件，在差分方程中加入对应的边界条件。
(5)迭代求解：根据差分方程和边界条件，进行迭代求解，更新网格上的数值。
(6)输出结果：根据需要，输出求解结果，可进行可视化显示或保存为文件。
3. Matlab中的实现在Matlab中，可以使用矩阵运算和循环结构来实现古典显式差分格式。
下面以一维热传导方程为例，介绍具体的实现步骤：(1)网格划分：通过定义离散网格的间隔dx，将求解区域划分为N 个小网格。
(2)时间步长确定：确定离散的时间步长Δt。
(3)差分方程建立：根据一维热传导方程，使用中心差分格式进行近似，得到差分方程。
(4)边界条件处理：根据问题的具体要求，确定边界条件。
(5)迭代求解：使用循环结构，依次求解每个时间步的网格节点的数值，并更新网格上的数值。
(6)输出结果：根据需要，可以将求解结果进行可视化显示或保存为文件。
在实际操作中，可以通过编写一个函数来实现古典显式差分格式的求解，函数的输入为问题的参数和边界条件，输出为求解的数值结果。
总结：古典显式差分格式是一种常用的数值方法，用于求解偏微分方程。
在Matlab中，可以通过矩阵运算和循环结构来实现古典显式差分格式的求解。
古典显式差分格式matlab
古典显式差分格式matlab古典显式差分格式（Classical Explicit Finite Difference Method）是一种常用的数值解法，用于求解偏微分方程（Partial Differential Equation, PDE）的近似解。
在本文中，我们将详细介绍古典显式差分格式的原理和步骤，并结合MATLAB代码进行实际应用。
1. 引言偏微分方程是一类重要的数学模型，用于描述许多科学和工程问题。
然而，求解偏微分方程的解析解通常较为困难甚至不可能，因此需要借助数值方法进行求解。
差分法是一种常见的数值方法，通过将坐标空间和时间空间离散化，将偏微分方程转化为代数方程，从而得到近似解。
古典显式差分格式是差分法中的一种常用方法，具有简单、易于实现的特点。
2. 原理古典显式差分格式的基本思想是对空间和时间进行离散化，将偏微分方程转化为差分方程。
具体而言，通过将空间和时间分割为离散的网格点，用离散点上的函数值近似代替偏微分方程中的导数项。
考虑一维热传导方程为例：∂u/∂t = α∂²u/∂x²其中u(x, t)为待求解的函数，α为常数。
我们需要求解该方程在给定初始条件和边界条件下的近似解。
为了离散化该方程，我们将空间和时间分别划分为N个小区间，得到N个网格点。
设x_i = iΔx，t_n = nΔt，其中Δx和Δt为空间和时间间隔。
引入差分近似，我们有以下显式差分格式：u_i^(n+1) = u_i^n + αΔt/Δx²* (u_(i-1)^n - 2u_i^n + u_(i+1)^n)其中u_i^n为时间步n和空间点i处的近似解，u_i^(n+1)为时间步n+1和空间点i处的近似解。
可以看出，当前时刻的解可以通过上一个时刻和邻近点的解计算得出。
3. 步骤为了使用古典显式差分格式求解偏微分方程，我们需要按照以下步骤进行操作：步骤1：确定空间和时间的离散化步长Δx和Δt，并计算网格节点数N。
差分方法基础
第二讲有限差分法基本原理一般的流体控制方程都是非线性的偏微分方程。
在绝大多数情况下，这些偏微分 方程无法得到精确解；而CFD 就是通过采用各种计算方法得到这些偏微分方程的数值 解，或称近似解。
当然这些近似解应该满足一定的精度。
目前，主要采用的CFD 方法是有限差分法和有限体积法。
本讲主要介绍有限差分法，它也是下一讲中的有限体积 法的基础［1］。
有限差分法求解流动控制方程的基本过程是：首先将求解区域划分为差分网格， 用有限个网格点代替连续的求解域，将待求解的流动变量(如密度、速度等)存储在 各网格点上，并将偏微分方程中的微分项用相应的差商代替，从而将偏微分方程转化 为代数形式的差分方程，得到含有离散点上的有限个未知变量的差分方程组。
求出该 差分方程组的解，也就得到了网格点上流动变量的数值解。
2.1差分和逼近误差由于通常数字计算机只能执行算术运算和逻辑运算，因此就需要一种用算术运算来处理函数微分运算的数值方法。
而有限差分法就是用离散网格点上的函数值来近似 导数的一种方法。
设有x 的解析函数y 二f(x)，从微分学知道函数y 对x 的导数为dy 、dx 分别是函数及自变量的微分，dy/dx 是函数对自变量的导数，又称微商。
相应地，上式中的 “、厶y 分别称为自变量及函数的差分，cy/^x 为函数对自变量的差 商。
在导数的定义中是以任意方式逼近于零的，因而.収是可正可负的。
在差分方 法中，“X 总是取某一小的正数。
这样一来，与微分对应的差分可以有三种形式：向前差分 二y = f (x 二x) _ f (x) 向后差分 y = f (x) _ f (x 中心差分11y 二 f (xx) _ f (xx)2 2上面谈的是一阶导数，对应的称为一阶差分。
对一阶差分再作一阶差分，就得到二阶差分，记为2y 。
以前向差分为例，有2A y =也(也y)=.'■： f (x _x ) - f (x) 1dx x —° x f(x：x) - f (x)Z x (2-1)=.f (x - x) - ■ :f(x) (2-2)-f (x 2 = x) 一f (x =x) I - If (x _x) - f (x) 1=f (x 2L X) -2 f (x =x) f (x) 依次类推，任何阶差分都可以由低一阶再作一阶差分得到。
中心差分格式 隐式 -回复
中心差分格式隐式-回复中心差分格式（Central Difference Scheme）是数值分析中一种常用的差分格式，常用于求解偏微分方程的数值解。
该格式的特点是在时间和空间上都采用中心差分的方式，能够提供较高的数值稳定性和精度。
在本文中，将详细介绍中心差分格式的基本原理、求解步骤以及应用范围。
中心差分格式的基本原理是利用有限差分法来近似求解偏微分方程。
偏微分方程通常描述了物理系统的变化过程，通过将连续的问题离散化为离散的数值问题，可以方便地使用计算机进行求解。
为了理解中心差分格式的原理，我们首先需要了解差分近似。
通过将时间和空间上的连续变量离散化，得到一系列有限的格点。
偏导数可以通过求取两个相邻格点之间的差商来近似计算。
在求解偏微分方程时，中心差分格式使用了中心差分来近似偏导数。
对于空间上的差分，中心差分将函数在当前格点和相邻格点上的函数值取平均，得到了近似的一阶导数。
对于时间上的差分，中心差分则利用了当前时刻和前一时刻的函数值来近似计算。
这种方式可以提供更高的数值稳定性，因为它能够考虑到当前时刻的变化以及前一时刻的信息。
中心差分格式的隐式形式是指使用未知量的未来值来近似计算当前值。
这样的形式可以提供更高的数值稳定性，但却需要通过一种迭代算法来求解。
下面我们将详细说明中心差分格式的求解步骤：1. 确定求解区域以及网格格点数量。
根据具体问题的要求，确定网格的大小和分布。
2. 初始化边界条件和初始条件。
根据实际情况，设置边界条件和初始条件。
3. 迭代求解。
通过迭代算法逐步更新未知量的值，直到达到收敛条件。
a) 在时间上进行迭代。
对于每个时间步长，更新时间点上未知量的值。
这里使用中心差分来近似计算。
b) 在空间上进行迭代。
对于每个空间点，更新其值。
这里同样使用中心差分来近似计算。
这个过程可以通过构建一个线性方程组来实现。
c) 判断是否达到收敛条件。
通常可以通过比较当前迭代结果和上一次迭代结果之间的差异来进行判断。
abaqus显示和隐式算法的差别
1.1. 弹性动力学有限元基本解法结构系统的通用运动学方程为：tR KU U C U M =++ （1） 求解该动力学振动响应主要有三类方法：（1）时域法（2）频域法（3）响应谱法 时域法又可分为：（1）直接积分法，（2）模态叠加法。
直接积分法又可分为中心差分法（显式），Wilson θ（隐式）法以及Newmark （隐式）法等。
本文介绍中心差分法（显式）与Newmark （隐式）法。
1 中心差分法（显式）假定0，1t ，2t ，…，n t 时刻的节点位移，速度与加速度均为已知，现求解)(t t t n ∆+时刻的 结构响应。
中心差分法对加速度，速度的导数采用中心差分代替，即为：)2(12t t t t t t U U U t U ∆+∆-+-∆= )(21t t t t t U U tU ∆-∆+-∆= (2) 将（2）式代入（1）式后整理得到tt t R U M ˆˆ=∆+ （3） 式（3）中C tM t M ∆+∆=211ˆ2 t t t t t U C tM t U M t K R R ∆-∆-∆-∆--=)211()2(ˆ22 分别称为有效质量矩阵，有效载荷矢量。
R ，M ，C ，K 为结构载荷，质量，阻尼，刚度矩阵。
求解线性方程组（3），即可获得t t ∆+时刻的节点位移向量t t U ∆+，将t t U ∆+代回几何方程与物理方程，可得t t ∆+时刻的单元应力和应变。
中心差分法在求解t t ∆+瞬时的位移t t U ∆+时，只需t t ∆+时刻以前的状态变量t U 和t t U ∆-，然后计算出有效质量矩阵M ˆ，有效载荷矢量tR ˆ，即可求出t t U ∆+，故称此解法为显式算法。
中心差分法，在开始计算时，需要仔细处理。
t =0时，要计算t U ∆，需要知道tU ∆-的值。
因此应该有一个起始技术，因而该算法不是自动起步的。
由于0U ，0U ，0U 是已知的，由t =0时的（2）式可知：02002U t U t U U t ∆+∆-=∆- 中心差分法中时间步长t ∆的选择涉及两个方面的约束：数值算法的稳定性和计算时间。
显式算法和隐式算法的介绍与比较
求解时间比较
使用显式方法，计算成本消耗与单元数量成正比，并且大致与最小 单元的尺寸成反比；分析中最消耗CPU的一项就是单元的处理。由于 CPU时间与积分的个数成正比，采用简化积分的单元便可极大地节省数 据存储量和运算次数，进而提高运算效率。但单点积分容易出现零能 模型且应力结果精确度与积分点相关。
（5）
U t t
Ut
Utt
[(1 2
)Ut Utt ]t2
（6）
将（5）（6）代入（1）中，得：
KˆUtt Rˆtt
其中：
Kˆ
1
t 2
M
t
C
K
（有效刚度矩阵）
Rˆtt
Rtt
1
M
[ t
2
Ut
1 t
Ut
(1 2
1)Ut ]
C[ t
U
t
(
1)Ut
( 2
1)tUt ]（有效载荷矢量）
可以看出，解当前 Utt，需要用到当前时刻 的 Rtt ，因此为隐式算法，需要用迭代实现。
特点：
当参数 0.5 ， 0.25(0.5 )2 时，Newmark法是无 条件稳定的，即 t 的大小不影响数值的稳定性。 此时时间步长的选择主要根据解的精度确定。
显式算法
总结
动态显式算法采用动力学方程的一些差分格式，不用直接求解切线 刚度，不需要进行平衡迭代，计算速度快，时间步长只要取的足够小， 一般不存在收敛性问题。因此需要的内存也比隐式算法要少。并且数值 计算过程可以很容易地进行并行计算，程序编制也相对简单。
LS-DYNA中中心差分法的步长控制
1D杆、梁单元：
PDE数值计算的有限差分法
PDE 数值计算的有限差分法《图像处理的PDE 方法》对给定的PDE 往往很难求其解析解，尤其是在实际问题中，这就需要求助于数值计算以获取该问题的近似解，常用的PDE 数值方法有有限差分法、有限元法和谱法等，其中，有限差分法应用得最为广泛。
因为待处理的图像通常已经是在二维空间中，按等采样而得到的离散化数字图像，这就自然构成了有限差分法所需要的等分网格（Grid ）。
1、有限差分格式有限差分的基本思想是：利用相距有限距离的两邻点的函数值的差与两点间距离的比值来近似函数对变量的偏导数。
例如，用向前差分来近似对时间的偏导数tu ∂∂，即n it nin iniuD tu u tu )(1++=∆-=∂∂对于空间中的一阶偏导数，除上面的向前差分外，还有向后差分、中心差分等，如下：向前差分：n ix nin i niuD tu u xu )(1++=∆-=∂∂后向差分：n ix ni n i niuD xu u xu )(1--=∆-=∂∂中心差分：n ix ni n i niuD xu u xu )0(112=∆-=∂∂-+根据泰勒展开式，有()+∆∂∂+∆∂∂+=∆+22221)()(x xu x xu x u x x u因此可得)()()(x O xx u x x u xu ∆+∆-∆+=∂∂说明向前差分和向后差分是一阶精度的。
同时，由于()+∆∂∂+∆∂∂+=∆+22221)()(x xu x x u x u x x u ()+∆∂∂+∆∂∂-=∆-22221)()(x xu x xu x u x x u可得())(2)()(2x O xx x u x x u xu ∆+∆∆--∆+=∂∂说明中心差分是二阶精度的。
当偏微分议程中含有二阶偏导数时，同样采用有限差分进行处理，先求出两个半点处的一阶偏导数中心差分，如下：x u u x u ni ni ni ∆-=⎥⎦⎤⎢⎣⎡∂∂++12/1，x u u x u ni n ini ∆-=⎥⎦⎤⎢⎣⎡∂∂--12/1然后再利用这两个一阶差分，作一次中心差分，得：()n ixx ni n i n i ni n i niuD x u u u xx u x u x u )0(2112/12/1222=∆+-=∆⎥⎦⎤⎢⎣⎡∂∂-⎥⎦⎤⎢⎣⎡∂∂=⎥⎦⎤⎢⎣⎡∂∂-+-+对于二阶偏导数yx u ∂∂∂2，同样采用类似的方法来处理，如下：x u u u u x u u x u j i j i j i j i j i j i j i ∆--+=∆-=⎥⎦⎤⎢⎣⎡∂∂+--++++-+++421,1,11,1,12/1,12/1,12/1,x u u u u x u u x u j i j i j i j i j i j i j i ∆--+=∆-=⎥⎦⎤⎢⎣⎡∂∂---+-+---+-42,11,1,11,12/1,12/1,12/1,其中()1,1,12/1,121++++++=j i ji j i u u u ()1,1,12/1,121+--+-+=j i ji j i u u u ()j i j i j i u u u ,11,12/1,121+-+-++= ()j i j i j i u uu ,11,12/1,121-----+=因此，yx u u u u yx u u u u u u u u yx u u u u yx u u u u yx u x u y x u j i j i j i j i ji j i j i j i j i j i j i j i ji j i j i j i j i j i j i j i j i j i nj i ∆∆--+=∆∆----+++=∆∆--+-∆∆--+=∆⎥⎦⎤⎢⎣⎡∂∂-⎥⎦⎤⎢⎣⎡∂∂=⎥⎦⎤⎢⎣⎡∂∂∂-++---+++-++-----+++---+-++--+++-+44441,11,11,11,1,11,11,1,1,11,11,1,1,11,1,11,11,1,11,1,12/1,2/1,,22、显式、隐式和半隐式方案以一维Burgers 方程xu ut u ∂∂=∂∂来说明几种PDE 的数值计算方案。
有限差分法
有限差分法有限差分法finite difference method微分方程和积分微分方程数值解的方法。
基本思想是把连续的定解区域用有限个离散点构成的网格来代替，这些离散点称作网格的节点；把连续定解区域上的连续变量的函数用在网格上定义的离散变量函数来近似；把原方程和定解条件中的微商用差商来近似，积分用积分和来近似，于是原微分方程和定解条件就近似地代之以代数方程组，即有限差分方程组，解此方程组就可以得到原问题在离散点上的近似解。
然后再利用插值方法便可以从离散解得到定解问题在整个区域上的近似解。
有限差分法的主要内容包括：如何根据问题的特点将定解区域作网格剖分；如何把原微分方程离散化为差分方程组以及如何解此代数方程组。
此外为了保证计算过程的可行和计算结果的正确，还需从理论上分析差分方程组的性态，包括解的唯一性、存在性和差分格式的相容性、收敛性和稳定性。
对于一个微分方程建立的各种差分格式，为了有实用意义，一个基本要求是它们能够任意逼近微分方程，这就是相容性要求。
另外，一个差分格式是否有用，最终要看差分方程的精确解能否任意逼近微分方程的解，这就是收敛性的概念。
此外，还有一个重要的概念必须考虑，即差分格式的稳定性。
因为差分格式的计算过程是逐层推进的，在计算第n＋1层的近似值时要用到第n层的近似值，直到与初始值有关。
前面各层若有舍入误差，必然影响到后面各层的值，如果误差的影响越来越大，以致差分格式的精确解的面貌完全被掩盖，这种格式是不稳定的，相反如果误差的传播是可以控制的，就认为格式是稳定的。
只有在这种情形，差分格式在实际计算中的近似解才可能任意逼近差分方程的精确解。
关于差分格式的构造一般有以下3种方法。
最常用的方法是数值微分法，比如用差商代替微商等。
另一方法叫积分插值法，因为在实际问题中得出的微分方程常常反映物理上的某种守恒原理，一般可以通过积分形式来表示。
此外还可以用待定系数法构造一些精度较高的差分格式。
高等油藏物理 第7章-差分格式的构造
du dt 0
u c
函数沿特征线不变
u(P) u(D ) u(B )
at x
u (B ) u ( A)
n 1
n
A P
特征线
用结点值形式写出
u
n 1 j
u
n j
a
t x
t
n
(u
n j
u
) j 1
D B
C
j j 1
n 2 j
u x x O x
n 3 j
4
一阶中心差商
u
n j 1
u
n j 1
2x
ux u
n
j
n
1 6
u x j
n 2
x O x
2
3
u x
j
n
u
n j 1
j 1
2x
O x
x
2
一阶偏导数的中心差商表达。它具有 阶的截断 误差，记为 ( x 2 ) 或者说距离有二阶精度。 R
n j
ux
j
n
x x
1 2 1 2
u x x
n 2 j
1 6 1 6
u x x O x
n 3 j
4
u j 1 u x j x , t n
n
u j u x
n
j
n
u x x
格式截断误差
R O t, x
j
j 1
（3）向后差分格式（左偏）
t n 1 n n n uj uj a u j u j 1 x u 0 f x j j