下载文档原格式
差分方法的原理和应用1. 原理介绍差分方法是一种数值计算方法，通过利用函数在某点附近的导数来近似计算函数的值。
差分方法主要基于以下两个原理：1.1 前向差分前向差分是通过计算函数在某点和其前面一个点的差值来近似计算函数的导数。
假设函数 f(x) 在点 x 处的导数为f’(x)，则前向差分的公式可以表示为：f'(x) ≈ (f(x+h) - f(x))/h其中，h 是一个小的正数，表示所选取的差分步长。
1.2 中心差分中心差分是通过计算函数在某点前后两个点的差值来近似计算函数的导数。
假设函数 f(x) 在点 x 处的导数为f’(x)，则中心差分的公式可以表示为：f'(x) ≈ (f(x+h) - f(x-h))/(2h)同样，h 是一个小的正数，表示所选取的差分步长。
2. 应用案例差分方法在许多科学和工程领域中都有广泛的应用。
以下列举了几个常见的应用案例：2.1 数值求导差分方法可以用于数值求导，即通过差分近似计算函数在某点处的导数。
通过选择合适的差分步长，可以获得足够高的精度。
数值求导在计算机图形学、数值分析等领域中被广泛使用。
2.2 数值积分差分方法还可以用于数值积分，即通过将函数离散化为一系列的差分点，然后计算这些差分点的和来近似计算函数的积分。
差分方法在求解常微分方程、偏微分方程等问题中也有重要的应用。
2.3 数据平滑差分方法可以用于数据平滑，即通过计算数据点之间的差分来减小数据的噪声。
通过选择合适的差分步长和平滑算法，可以过滤掉数据中的噪声，并提取出数据的趋势。
2.4 图像处理差分方法在图像处理中也有广泛的应用。
例如，图像边缘检测算法就是基于差分方法的。
通过计算图像中像素之间的差分，可以检测出图像中的边缘。
2.5 数值优化差分方法还可以用于数值优化，即通过利用函数在某点附近的差分信息来搜索函数的最优解。
差分方法在机器学习、优化算法中有重要的应用。
3. 总结差分方法是一种常见的数值计算方法，通过利用函数在某点附近的导数来近似计算函数的值。
中心差分法的基本理论与程序设计1程序设计的目的与意义该程序通过用C语言（部分C++语言）编写了有限元中用于求解动力学问题的中心差分法，巩固和掌握了中心差分法的基本概念，提高了实际动手能力，并通过实际编程实现了中心差分法在求解某些动力学问题中的运用，加深了对该方法的理解和掌握。
2程序功能及特点该程序采用C语言（部分C++语言）实现了用于求解动力学问题的中心差分法，可以求解得到运动方程的解答，包括位移，速度和加速度。
计算简便且在算法稳定的条件下，精度较高。
3中心差分法的基本理论在动力学问题中，系统的有限元求解方程（运动方程）如下所示：阵和结点载荷向量，并分别由各自的单元矩阵和向量集成。
与静力学分析相比，在动力分析中，由于惯性力和阻尼力出现在平衡方程中，因此引入了质量矩阵和阻尼矩阵，最后得到的求解方程不是代数方程组，而是常微分方程组。
常微分方程的求解方法可以分为两类，即直接积分法和振型叠加法。
中心差分法属于直接积分法，其对运动方程不进行方程形式的变换而直接进行逐步数值积分。
通常的直接积分是基于两个概念，中心差分法的基本思路是用有限差分代替位移对时间的求导，将运动方程中的速度和加速度用位移的某种组合表示，然后将常微分方程组的求解问题转换为代数方程组的求解问题，并假设在每个小的时间间隔内满足运动方程，则可以求得每个时间间隔的递推公式，进而求得整个时程的反应。
在中心差分法中，加速度和速度可以用位移表示，即：(2t =∆(2t=∆而得到。
为此将加速度和速度的表达式代入上式中，即可得到中心差分法的递推公式：个离散时间点的解的递推公式，这种数值积分方法又称为逐步积分法。
需要指出和速度的表达式可知：2a +中心差分法避免了矩阵求逆的运算，是显式算法，且其为条件稳定算法，利型载荷引起的波传播问题的求解，而对于结构动力学问题则不太合适。
4 中心差分法的有限元计算格式利用中心差分法逐步求解运动方程的算法步骤如下所示： 1. 初始计算(1)(2)(3)ω=；(4)(5)(6)2.(1)(2)(3)5 程序设计5.1 程序流程图1 程序流程图各子程序主要功能为：ArrayLU ：LU Inverse ：求矩阵的转置矩阵； ArrayMVector ：矩阵和向量的乘法； LUSolve5.2 输入数据及变量说明5.2.1 输入数据该程序的原始输入数据应包括三个部分： (1)(2)(3) 确定时间步长，其中为了保证该算法的稳定性，需要满足ω=5.2.2 变量说明该程序的各个变量含义如下： (1) num ，timeStep ，dtnum ——矩阵维度； timeStep ——时间步数； dt ——时间步长；(2) M ，C ，K ，X ，V ，A ，P ，MM ，PT ，c0，c1，c2，c3M ——质量矩阵； C ——阻尼矩阵； K ——刚度矩阵； X ——位移矩阵； V ——速度矩阵；A ——加速度矩阵； P ——载荷向量； MM ——有效质量矩阵； PT —— c0,c1,c2,c3——积分常数；6算例6.1问题描述应用本程序计算一个三自由度系统，它的运动方程是：⎣当时，利用公式，可以计算得到6.2理论计算6.2.1中心差分法（理论解）(1)则起步条件为⎢⎥⎣⎦对于每一时间步长，需先计算有效载荷：由上式得到的每一时间步长的位移结果如表1所示：表10.00 0.00 0.00 0.03 0.13 0.36 0.79 1.46 2.373.42 0.00 0.03 0.19 0.58 1.26 2.24 3.434.695.846.770.401.482.974.525.826.717.227.517.858.45(2)按照相同的步骤，所得结果如下：再计算下去，位移将继续增大，这是不稳定的典型表现。
两点中心差分公式推导
两点中心差分公式是数值分析中的一种重要公式，用于计算两个连续变量之间的差分。
该公式可以表示为：
Δy = (y2 - y1) / (x2 - x1)
其中，Δy表示两点之间的纵向差分，y1和y2分别表示两点中的纵坐标值，x1和x2分别表示两点中的横坐标值。
下面是两点中心差分公式的推导过程：
1. 首先，设点1的坐标为(x1, y1)，点2的坐标为(x2, y2)。
2. 计算两点之间的纵向差分Δy = y2 - y1。
3. 计算两点之间的横坐标差分Δx = x2 - x1。
4. 将Δy和Δx代入公式Δy = (y2 - y1) / (x2 - x1)，得到两点中心差分公式。
另外，两点中心差分公式也可以用于计算平均变化率。
当Δx不为零时，两点中心差分公式可以表示为：
Δy/Δx = (y2 - y1) / (x2 - x1)
这表示了变量y在两点间的平均变化率。
注意，当Δx为零时，两点中心差分公式不适用，因为分母不能为零。
在这种情况下，可以采用其他方法（如前后差分法）来计算变化率。
中心差分法计算程序编程姓名：张泽伟 学号： 电话：一、中心差分法程序原理说明1。
1 中心差分法思路中心差分法的基本思路：是将运动方程中的速度向量和加速度向量用位移的某种组合来表示，将微分方程组的求解问题转化为代数方程组的求解问题，并在时间区间内求得每个微小时间区间的递推公式,进而求得整个时程的反应。
1。
2 中心差分法原理中心差分法只在相隔t ∆一些离散的时间区间内满足运动方程，其基于有限差分代替位移对时间的求导(即速度和加速度)，如果采用等时间步长，t t i ∆=∆,则速度与加速度的中心差分近似为：t u u u i i ∆-=-+•211 (a ） 2112t u u u u i i i ∆+-=-+•• （b ） 而离散时间点的运动为)(),(),(i i i i i i t u u t u u t u u ••••••=== ( =i 0，1,2,3,……）由体系运动方程为：0)()()(=++•••t ku t u c t u m i (c) 将速度和加速度的差分近似公式（a ）和式(b)代入式(c ）可以得到i t 时刻的运动方程：02211211=+∆-+∆+--+-+i i i i i i ku t u u c t u u u m （d ）在（d ）式中，假设i u 和1-i u 是已知的，即在i t 及i t 以前时刻的运动已知，则可以把已知项移到方程的右边，整理得到：12212)2()2()2(-+∆-∆-∆--=∆+∆i i i u t c t m u t m k u t c t m (e)由式（e ）就可以根据i t 及i t 以前时刻的运动，求得1+i t 时刻的运动，如果需要可以用式(a ）和式（b)求得体系的速度和加速度。
1.3 初始条件转化假设给定的初始条件为),0(),0(00••==u u u u （g ）由式（g ）确定1-u .在零时刻速度和加速度的中心差分公式为：t u u u ∆-=-•2110 (h)` 210102t u u u u ∆+-=-•• （i)将式(i ）消去1u 得:020012•••-∆+∆-=u t u t u u （j ）而零时刻的加速度值0••u 可以用t ＝0时的运动方程0000=++•••ku u c u m确定 即 )(1000ku u c m u --=••• （k ） 这样就可以根据初始条件00,•u u 和初始荷载0P ，就可以根据上式确定1-u 的值。
中心差分法的基本理论与程序设计
中心差分法（Central Difference Method）是一种常用的数值计算方法，用于近似求解函数的导数。
该方法基于导数的定义，通过两个函数值的差分来近似求解导数的值。
本文将介绍中心差分法的基本理论和程序设计过程。
一、中心差分法的基本理论
f'(x)≈(f(x+h)-f(x-h))/(2h)
其中，h是步长，决定了计算函数值的邻近点的间距。
二、中心差分法的程序设计
下面我们将介绍中心差分法的程序设计过程，包括步骤和代码实现。
步骤：
1.定义函数f(x)，表示需要求解导数的函数。
2.定义步长h，决定计算函数值的邻近点的间距。
3.根据中心差分法的公式，计算导数的近似值：
f'(x)≈(f(x+h)-f(x-h))/(2h)
4.返回导数的近似值。
代码实现：
```
def central_difference(f, x, h):
"""
计算函数f在点x处的导数值
:param f: 函数f(x)
:param x: 需要求导数的点
:param h: 步长
:return: 导数的近似值
"""
df = (f(x + h) - f(x - h)) / (2 * h)
return df
```
以上就是中心差分法的基本理论和程序设计过程。
中心差分法是一种简单有效的数值计算方法，可以用于求解导数的近似值。
在实际应用中，可以根据需要选择合适的步长，并通过增加点的数量来提高计算结果的精度。
中心差分法计算单自由度体系的动力反应北京工业大学结构工程组员:胡建华 S201204111马 恒 S201204112 陈相家S201204083 张力嘉S201204022 0前言时域逐步积分法是数值分析方法，它只假设结构本构关系在一个微小的时间步距内是线性的。
时域逐步积分法是结构动力问题中一个得到广泛研究的课题，并在结构动力反应计算中得到广泛应用。
由于引进的假设条件不同，可以有各种不同的方法，比如中心差分法，线性加速度法，平均常加速度法，Wilson-θ法等，其中中心差分法精度最高。
在本文中，通过编制中心差分法计算单自由度体系的动力反应的程序来了解其应用及稳定性。
1中心差分法原理中心差分法的基本思路：是将运动方程中的速度向量和加速度向量用位移的某种组合来表示，将微分方程组的求解问题转化为代数方程组的求解问题，并在时间区间内求得每个微小时间区间的递推公式，进而求得整个时程的反应。
中心差分法只在相隔t ∆一些离散的时间区间内满足运动方程，其基于有限差分代替位移对时间的求导（即速度和加速度），如果采用等时间步长，t t i ∆=∆，则速度与加速度的中心差分近似为： tu u u i i ∆-=-+•211 (a)2112tu u u u i i i ∆+-=-+•• (b) 而离散时间点的运动为)(),(),(i i i i i i t u u t u u t u u ••••••=== ( =i 0,1,2,3,……) 由体系运动方程为：0)()()(=++•••t ku t u c t u m i (c)将速度和加速度的差分近似公式（a ）和式（b ）代入式（c ）可以得到i t 时刻的运动方程：02211211=+∆-+∆+--+-+i i i i i i ku t u u c t u u u m（d ）在（d ）式中，假设i u 和1-i u 是已知的，即在i t 及i t 以前时刻的运动已知，则可以把已知项移到方程的右边，整理得到：12212)2()2()2(-+∆-∆-∆--=∆+∆i i i u t c tm u t m k u t c t m (e) 由式（e ）就可以根据i t 及i t 以前时刻的运动，求得1+i t 时刻的运动，如果需要可以用式（a ）和式（b ）求得体系的速度和加速度。
中心差分法计算程序编程具体而言，中心差分法首先选择一个较小的步长h，然后计算对应点的函数值。
接着，利用这两个较近的点对函数在该点处的导数进行逼近，即利用函数值的差分来计算导数的差分。
通过这种方式，可以得到一个近似的导数值。
以下是一个用Python编写的中心差分法计算程序示例：```pythonimport numpy as npdef center_diff(f, x, h):""""""return (f(x + h) - f(x - h)) / (2 * h)def func(x):"""待求导函数"""return np.sin(x)def main(:x=1.0#需要求导的点h=0.1#步长derivative = center_diff(func, x, h)print("导数的近似值：", derivative)if __name__ == "__main__":main```在这个示例程序中，`center_diff` 函数接受一个函数 `f`、一个点`x` 和一个步长 `h` 作为参数，计算函数 `f` 在点 `x` 处的导数。
`func` 函数定义了要求导的函数，这里以求 sin 函数的导数为例。
在 `main` 函数中，定义了需要求导的点 `x` 和步长 `h`，然后调用 `center_diff` 函数计算导数的近似值，并将结果打印出来。
使用该程序，可以将中心差分法用于各种函数的求导计算。
只需要将待求导的函数定义为 `func` 函数，并根据需要调整求导点 `x` 和步长`h`。
需要注意的是，步长`h`的选择要适中。
如果`h`太大，差分近似可能不准确；如果`h`太小，计算时会引入较大的数值误差。
通常，可以尝试不同的`h`值，选择一个能够在精度和计算效率之间取得平衡的值。
差分法基本原理范文差分法是一种数值求解微分方程的方法，它将微分方程转化为离散形式，通过离散点上的函数值之差来近似计算导数，从而得到微分方程的数值解。
差分法的基本原理可以总结为以下几个步骤：1.网格划分：将求解区域划分为若干有限个小区域，每个小区域称为一个网格单元。
通常采用均匀网格划分，将区域划分为有限个等距的小区域。
2.网格节点：在每个网格单元的边界上选择一个或多个节点，节点是离散区域内的点。
节点数量取决于所选择的差分格式。
通常要求节点密度足够高，以确保数值解的精度。
3.差分逼近：使用差分公式对微分方程中的导数进行近似。
差分公式的选择是差分法的核心。
常见的差分公式包括：中心差分、向前差分、向后差分等等。
不同的差分公式对应着不同的差分格式，如前向差分格式、后向差分格式、中心差分格式等等。
4.离散化方程：根据差分逼近的结果，将微分方程中的导数用离散点上的函数值之差来代替，得到离散的差分方程。
离散化的过程将微分方程转化为代数方程组，可以通过求解代数方程组来得到数值解。
5.边界条件：确定边界条件，在差分方程中将其加入到方程组中。
边界条件通常是指在求解区域边界上的已知函数值或导数值。
6.求解代数方程组：根据离散化方程和边界条件，得到一个代数方程组。
通过数值方法，如高斯消元法、迭代法等，求解得到方程组的数值解。
7.结果输出和误差分析：根据求解得到的数值解，可以进行结果输出、误差分析和收敛性研究。
通常需要对数值解进行采样，与解析解进行比较，以评估差分法的精度和稳定性。
差分法的优点包括：简单易用，计算效率高，适用于各种类型的微分方程，比如常微分方程、偏微分方程及椭圆、抛物、双曲型方程等等。
然而，差分法也存在一些限制，如数值方法的稳定性与精度受节点密度、步长选择和差分格式的影响。
因此，在实际应用中，需要根据具体问题的特点和求解要求选择合适的差分格式和参数，以获得满足要求的数值解。
中心差分法的基本理论与程序设计1程序设计的目的与意义该程序通过用C语言（部分C++语言）编写了有限元中用于求解动力学问题的中心差分法，巩固和掌握了中心差分法的基本概念，提高了实际动手能力，并通过实际编程实现了中心差分法在求解某些动力学问题中的运用，加深了对该方法的理解和掌握。
2程序功能及特点该程序采用C语言（部分C++语言）实现了用于求解动力学问题的中心差分法，可以求解得到运动方程的解答，包括位移，速度和加速度。
计算简便且在算法稳定的条件下，精度较高。
3中心差分法的基本理论在动力学问题中，系统的有限元求解方程（运动方程）如下所示：()()()()Ma t Ca t Ka t Q t++=式中，()a t分别是系统的结点加速度向a t是系统结点位移向量，()a t和()量和结点速度向量，,,M C K和()Q t分别是系统的质量矩阵、阻尼矩阵、刚度矩阵和结点载荷向量，并分别由各自的单元矩阵和向量集成。
与静力学分析相比，在动力分析中，由于惯性力和阻尼力出现在平衡方程中，因此引入了质量矩阵和阻尼矩阵，最后得到的求解方程不是代数方程组，而是常微分方程组。
常微分方程的求解方法可以分为两类，即直接积分法和振型叠加法。
中心差分法属于直接积分法，其对运动方程不进行方程形式的变换而直接进行逐步数值积分。
通常的直接积分是基于两个概念，一是将在求解域0t T内的任何时刻t都应满足运动方程的要求，代之仅在一定条件下近似地满足运动方程，例如可以仅在相隔t∆的离散的时间点满足运动方程；二是在一定数目的t∆区域内，假设位移a、速度a、加速度a的函数形式。
中心差分法的基本思路是用有限差分代替位移对时间的求导，将运动方程中的速度和加速度用位移的某种组合表示，然后将常微分方程组的求解问题转换为代数方程组的求解问题，并假设在每个小的时间间隔内满足运动方程，则可以求得每个时间间隔的递推公式，进而求得整个时程的反应。
在中心差分法中，加速度和速度可以用位移表示，即：()212t t t t t a a a a t -∆+∆=-+∆ ()12t t t t a a a t-∆+∆=-+∆时间t t ∆+的位移解答t t a +∆，可由时间t 的运动方程应得到满足，即由下式：t t t t Ma Ca Ka Q ++=而得到。
中心差分法计算程序编程姓名：张泽伟 学号： 电话：一、中心差分法程序原理说明 1.1 中心差分法思路中心差分法的基本思路：是将运动方程中的速度向量和加速度向量用位移的某种组合来表示，将微分方程组的求解问题转化为代数方程组的求解问题，并在时间区间内求得每个微小时间区间的递推公式，进而求得整个时程的反应。
1.2 中心差分法原理中心差分法只在相隔t ∆一些离散的时间区间内满足运动方程，其基于有限差分代替位移对时间的求导（即速度和加速度），如果采用等时间步长，t t i ∆=∆，则速度与加速度的中心差分近似为：t u u u i i ∆-=-+•211 (a) 2112t u u u u i i i ∆+-=-+•• (b)而离散时间点的运动为)(),(),(i i i i i i t u u t u u t u u ••••••=== ( =i 0,1,2,3,……)由体系运动方程为：)()()(=++•••t ku t u c t u m i (c)将速度和加速度的差分近似公式（a ）和式（b ）代入式（c ）可以得到i t时刻的运动方程：2211211=+∆-+∆+--+-+i i i i i i ku t u u c t u u u m（d ）在（d ）式中，假设i u 和1-i u 是已知的，即在i t 及i t以前时刻的运动已知，则可以把已知项移到方程的右边，整理得到：12212)2()2()2(-+∆-∆-∆--=∆+∆i i i u t c t m u t m k u t c t m (e)由式（e ）就可以根据i t 及i t 以前时刻的运动，求得1+i t时刻的运动，如果需要可以用式（a ）和式（b ）求得体系的速度和加速度。
1.3 初始条件转化假设给定的初始条件为),0(),0(00••==u u u u （g ） 由式（g ）确定1-u 。
在零时刻速度和加速度的中心差分公式为：t u u u ∆-=-•2110 （h ） `210102t u u u u ∆+-=-•• （i ）将式（i ）消去1u 得：020012•••-∆+∆-=u t u t u u （j ）而零时刻的加速度值0••u 可以用t ＝0时的运动方程 0000=++•••ku u c u m确定即 )(1000ku u c m u --=••• （k ）这样就可以根据初始条件00,•u u 和初始荷载P ，就可以根据上式确定1-u 的值。
中心差分格式隐式-回复中心差分格式是一种数值计算的方法，用于近似求解偏微分方程。
隐式中心差分格式是其中一种扩展形式，将迭代过程中的中心差分格式变为隐式求解，具有更好的稳定性和精确性。
本文将逐步介绍中心差分格式和隐式中心差分格式的原理、应用和算法实现过程。
一、中心差分格式的原理中心差分格式是一种通过近似导数来表示微分项的方法。
它的基本思想是以离散的方式计算函数在某一点的导数，通过将函数在该点展开为泰勒级数的形式，利用离散点的函数值进行近似计算。
中心差分表示的导数形式为：f'(x) ≈(f(x+Δx) - f(x-Δx)) / (2Δx)其中，f(x+Δx)和f(x-Δx)表示在x+Δx和x-Δx处的函数值，Δx表示离散点之间的间距。
二、隐式中心差分格式的原理隐式中心差分格式是中心差分格式的一种扩展形式，它的特点是在迭代过程中使用隐式求解。
隐式求解是指在迭代过程中，未知量出现在方程的两侧，需要通过求解线性方程组来求得。
隐式中心差分格式的迭代公式为：un+1 = un + λ(u(n+1) - 2un + u(n-1))其中，un表示时间步t=n时的解，un+1表示时间步t=n+1时的解，u(n-1)表示时间步t=n-1时的解，λ表示时间步长的比例系数。
三、隐式中心差分格式的应用隐式中心差分格式广泛应用于偏微分方程的数值解法中，特别是对于具有稳定性要求的方程。
由于其隐式求解的特点，可以在一定程度上提高数值解的稳定性，并且能够处理非线性项的情况。
隐式中心差分格式的应用领域包括热传导方程、扩散方程、波动方程等。
在这些方程中，隐式中心差分格式能够提供更为准确和稳定的数值解，并且具有更大的适用范围。
四、隐式中心差分格式的算法实现过程隐式中心差分格式的算法实现过程主要包括以下几个步骤：1. 网格的划分：将求解区域划分为离散的网格点，确定离散点之间的间距Δx和时间步长Δt。
2. 边界条件的处理：根据具体问题的边界条件，对应的确定边界点的数值。
中心差分格式隐式-回复什么是中心差分格式？如何使用中心差分格式进行数值计算？中心差分格式是一种常用的数值计算方法，用于求解偏微分方程的数值解。
它的基本思想是将连续的物理空间离散化为离散的网格点，并以离散方式逼近连续的微分算子。
中心差分格式通过使用当前节点和相邻节点的函数值来估计函数的导数，以实现数值计算。
中心差分格式的基本公式如下：\frac{{df}}{{dx}} \approx \frac{{f(x + \Delta x) - f(x - \Delta x)}}{{2\cdot \Delta x}}其中，\Delta x是网格间距，表示相邻节点之间的距离。
根据该公式，我们可以通过给定函数在相邻网格点的函数值来估计该点处的导数。
这种方法通常被称为“二阶中心差分”。
在使用中心差分格式求解偏微分方程时，我们通常将方程中的导数替换为中心差分公式，并将方程离散为离散点上的代数方程组。
通过解这个方程组，我们可以得到方程的数值近似解。
具体来说，使用中心差分格式进行数值计算的步骤如下：1. 确定计算区域和网格：首先需要确定计算区域的边界条件，并将计算区域划分为离散的网格点。
通常情况下，我们会选择等距离的网格点，使得每个网格点之间的间距相等。
2. 确定离散化的时间步长：在时间方向上，我们也需要进行离散化。
需要确定时间步长\Delta t，表示相邻时间层之间的间隔。
通常情况下，我们会选择足够小的时间步长，以保证数值解的稳定性和精度。
3. 计算方程的右侧项：将偏微分方程中的导数替换为中心差分公式，可以得到离散点上方程的代数形式。
这样，我们就可以计算方程的右侧项，包括边界条件和源项。
4. 构建代数方程组：根据离散后的方程形式，我们可以将方程表示为一个代数方程组。
通常情况下，我们会使用矩阵表示方程组，并利用矩阵运算求解代数方程组。
5. 求解代数方程组：通过求解代数方程组，我们可以得到方程的数值解。
通常情况下，我们会使用迭代方法（如雅可比迭代、高斯-赛德尔迭代等）或直接求解方法（如LU分解、共轭梯度法等）来求解代数方程组。
中心差分法的基本理论与程序设计1程序设计的目的与意义该程序通过用C语言（部分C++语言）编写了有限元中用于求解动力学问题的中心差分法，巩固和掌握了中心差分法的基本概念，提高了实际动手能力，并通过实际编程实现了中心差分法在求解某些动力学问题中的运用，加深了对该方法的理解和掌握。
2程序功能及特点该程序采用C语言（部分C++语言）实现了用于求解动力学问题的中心差分法，可以求解得到运动方程的解答，包括位移，速度和加速度。
计算简便且在算法稳定的条件下，精度较高。
3中心差分法的基本理论在动力学问题中，系统的有限元求解方程（运动方程）如下所示：()()()()Ma t Ca t Ka t Q t++=式中，()a t分别是系统的结点加速度向a t是系统结点位移向量，()a t和()量和结点速度向量，,,M C K和()Q t分别是系统的质量矩阵、阻尼矩阵、刚度矩阵和结点载荷向量，并分别由各自的单元矩阵和向量集成。
与静力学分析相比，在动力分析中，由于惯性力和阻尼力出现在平衡方程中，因此引入了质量矩阵和阻尼矩阵，最后得到的求解方程不是代数方程组，而是常微分方程组。
常微分方程的求解方法可以分为两类，即直接积分法和振型叠加法。
中心差分法属于直接积分法，其对运动方程不进行方程形式的变换而直接进行逐步数值积分。
通常的直接积分是基于两个概念，一是将在求解域0t T内的任何时刻t都应满足运动方程的要求，代之仅在一定条件下近似地满足运动方程，例如可以仅在相隔t∆的离散的时间点满足运动方程；二是在一定数目的t∆区域内，假设位移a、速度a、加速度a的函数形式。
中心差分法的基本思路是用有限差分代替位移对时间的求导，将运动方程中的速度和加速度用位移的某种组合表示，然后将常微分方程组的求解问题转换为代数方程组的求解问题，并假设在每个小的时间间隔内满足运动方程，则可以求得每个时间间隔的递推公式，进而求得整个时程的反应。
在中心差分法中，加速度和速度可以用位移表示，即：()212t t t t t a a a a t -∆+∆=-+∆ ()12t t t t a a a t-∆+∆=-+∆时间t t ∆+的位移解答t t a +∆，可由时间t 的运动方程应得到满足，即由下式：t t t t Ma Ca Ka Q ++=而得到。
为此将加速度和速度的表达式代入上式中，即可得到中心差分法的递推公式：2221121122t t t tt t M C a Q K M a M C a t t t t t +∆-∆⎛⎫⎛⎫⎛⎫+=---- ⎪ ⎪ ⎪∆∆∆∆∆⎝⎭⎝⎭⎝⎭若已经求得t a 和t t a -∆，则从上式可以进一步解出t t a +∆。
所以上式是求解各个离散时间点的解的递推公式，这种数值积分方法又称为逐步积分法。
需要指出的是，此算法有一个起步问题。
因为当0t =时，为了计算t a ∆，除了知道初始条件已知的0a ，还需要知道t a -∆，所以必须用一专门的起步方法。
根据以上加速度和速度的表达式可知：20002tt a a ta a -∆∆=-∆+其中0a 和0a 可以从给定的初始条件中得到，而0a 则可以利用0t =时的运动方程0000Ma Ca Ka Q ++=得到，即：()10000a M Q Ca Ka -=--中心差分法避免了矩阵求逆的运算，是显式算法，且其为条件稳定算法，利用它求解具体问题时，时间步长t ∆必须小于由该问题求解方程性质所决定的某个临界值cr t ∆，否则算法将是不稳定的。
中心差分法比较适用于由冲击、爆炸类型载荷引起的波传播问题的求解，而对于结构动力学问题则不太合适。
4 中心差分法的有限元计算格式利用中心差分法逐步求解运动方程的算法步骤如下所示： 1. 初始计算(1) 形成刚度矩阵K 、质量矩阵M 和阻尼矩阵C ； (2) 给定0a ，0a 和0a ； (3) 选择时间步长t ∆，2ncr nT tt ωπ∆∆==，并计算积分常数021c t =∆，112c t=∆，202c c =， 321c c =；(4) 计算0030t a a ta c a -∆=-∆+；(5) 形成有效质量矩阵01ˆM c M c C =+； (6) 三角分解ˆM：ˆM LU =。
2. 对于每一时间步长（0,,2t t t =∆∆⋅⋅⋅） (1) 计算时间t 的有效载荷()()201ˆt t t t tQ Q K c M a c M c C a -∆=---- (2) 求解时间t t +∆的位移ˆt t tLUa Q +∆= (3) 如果需要，计算时间t 的加速度和速度()02t t t t t t a c a a a -∆+∆=-+()1t t t t t a c a a -∆+∆=-+5程序设计5.1程序流程图1 程序流程图各子程序主要功能为：=；ArrayLU：LU三角分解A LUAInverse：求矩阵的转置矩阵；ArrayMVector：矩阵和向量的乘法；=。
LUSolve：求解方程LUX P5.2 输入数据及变量说明5.2.1 输入数据该程序的原始输入数据应包括三个部分： (1) 刚度矩阵K ，质量矩阵M 和阻尼矩阵C ；(2) 初始条件：时间0t =时刻的位移0a ，速度0a ，加速度0a ；(3) 确定时间步长t ∆，其中为了保证该算法的稳定性，需要满足2ncr nT tt ωπ∆∆==。
5.2.2 变量说明该程序的各个变量含义如下： (1) num ，timeStep ，dtnum ——矩阵维度； timeStep ——时间步数； dt ——时间步长；(2) M ，C ，K ，X ，V ，A ，P ，MM ，PT ，c0，c1，c2，c3M ——质量矩阵； C ——阻尼矩阵； K ——刚度矩阵； X ——位移矩阵； V ——速度矩阵； A ——加速度矩阵； P ——载荷向量； MM ——有效质量矩阵； PT ——时间t 时刻的有效载荷； c0,c1,c2,c3——积分常数；6 算例6.1 问题描述应用本程序计算一个三自由度系统，它的运动方程是：100210003014200010226a a -⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥+--=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥-⎣⎦⎣⎦⎣⎦初始条件：当0t =时，00a =，00a =。
已知此系统的固有频率为：1ω=2ω=3ω=110.89T =，2 4.444T =，3 3.628T =。
当0t =时，利用公式()10000a M Q Ca Ka -=--，可以计算得到[]0006Ta =；时间步长分别取3100.363t T ∆==和3518.14t T ∆==进行计算。
6.2 理论计算6.2.1 中心差分法（理论解）(1) 时间步长3100.363t T ∆==当3100.363t T ∆==时，其积分常数为：0217.589c t ==∆ 111.3772c t==∆ 20215.178c c == 321 6.5882c c e ==- 则起步条件为0030000000.36300.0659000060.3953ta a ta c a -∆⎡⎤⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥=-∆+=-+=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦⎣⎦有效质量矩阵ˆM为：011000007.5900ˆ7.59030 1.38000022.770001000007.59M c M c C ⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥=+=+=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦对于每一时间步长，需先计算有效载荷：()()201ˆ013.18107.59000141.532022.77060213.18007.59t t t t tt t tQ Q K c M a c M c C a a a -∆-∆=----⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥=+-⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦再从下列方程计算t t +∆时间的位移t t a +∆7.5900ˆ022.770007.59t t a Q +∆⎡⎤⎢⎥=⎢⎥⎢⎥⎣⎦由上式得到的每一时间步长的位移结果如表1所示：表1 3100.363t T ∆==理论解时间t ∆ 2t ∆ 3t ∆ 4t ∆ 5t ∆ 6t ∆ 7t ∆ 8t ∆ 9t ∆ 10t ∆ 1a 0.00 0.00 0.00 0.03 0.13 0.36 0.79 1.46 2.37 3.42 2a 0.00 0.03 0.19 0.58 1.26 2.24 3.43 4.69 5.84 6.77 3a0.401.482.974.525.826.717.227.517.858.45(2) 时间步长3518.14t T ∆==按照相同的步骤，所得结果如下：2009.8710t a ∆⎡⎤⎢⎥=⎢⎥⎢⎥⨯⎣⎦ 52502.07106.4610t a ∆⎡⎤⎢⎥=⨯⎢⎥⎢⎥⨯⎣⎦ 78387.13102.36105.6610t a ∆⎡⎤⨯⎢⎥=⨯⎢⎥⎢⎥⨯⎣⎦再计算下去，位移将继续增大，这是不稳定的典型表现。
其原因是在条件稳定的中心差分方法中采用了远大于()cr t T π∆=的时间步长()3518.14t T ∆==，所以不可能得到有意义的结果。
6.2.2 振型叠加法（精确解）采用振型叠加法对上述问题进行计算，可以得到该问题的精确解。
首先应求解的广义特征值问题为：2210100142030022001φωφ-⎡⎤⎡⎤⎢⎥⎢⎥--=⎢⎥⎢⎥⎢⎥⎢⎥-⎣⎦⎣⎦按照一般的线性代数方法可以得到上式的解答为：2113ω= 15123Tφ⎡⎤=⎢⎥⎣⎦222ω= []1201Tφ=-233ω= []1112Tφ=--利用上式，可以将原问题转换为以123,,φφφ为基向量的3个互不耦合的运动方程，即：()()1119103x t x t +=()()22265x t x t += ()()33332x t x t +=-原系统的初始条件是00a =，00a =，经转换后为：00i t x == 00i t x == ()1,2,3i =利用无阻尼情形的杜哈美积分公式，可以得到上述方程的解答为：())1101cos 27x t ⎡⎤=-⎣⎦ ())231cos 5x t ⎡⎤=-⎣⎦())311cos 2x t ⎡⎤=--⎣⎦最后利用振型叠加得到系统的位移为：())))101cos 27121353011cos 521211cos 2a t ⎡⎤⎡⎤-⎢⎥⎣⎦--⎡⎤⎢⎥⎢⎥⎢⎥⎡⎤=-⎢⎥⎣⎦⎢⎥⎢⎥-⎢⎥⎣⎦⎡⎤⎢⎥--⎣⎦⎢⎥⎣⎦根据上式计算得到的每一时间步长的位移值是系统响应的精确解，如表2 和表3所示。
(1) 3100.363t T ∆==的精确位移值。
表2 3100.363t T ∆==精确解时间t ∆ 2t ∆ 3t ∆ 4t ∆ 5t ∆ 6t ∆ 7t ∆ 8t ∆ 9t ∆ 10t ∆ 1a 0.00 0.00 0.01 0.04 0.14 0.37 0.79 1.45 2.32 3.34 2a 0.00 0.04 0.21 0.59 1.25 2.21 3.38 4.63 5.80 6.76 3a0.391.452.924.485.816.747.297.607.928.46(2) 3518.14t T ∆==的精确位移值。
