下载文档原格式
中心差分法的基本理论与程序设计
中心差分法（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
```
以上就是中心差分法的基本理论和程序设计过程。
中心差分法是一种简单有效的数值计算方法，可以用于求解导数的近似值。
在实际应用中，可以根据需要选择合适的步长，并通过增加点的数量来提高计算结果的精度。
中心差分动力学方程 matlab（实用版）目录1.引言2.中心差分的概念和应用3.动力学方程的概述4.MATLAB 在中心差分和动力学方程中的应用5.结论正文一、引言在科学研究和工程技术中，数学模型和计算方法的结合对于解决实际问题具有重要意义。
中心差分和动力学方程是数学模型和计算方法中的两个重要概念，它们在许多领域中都有着广泛的应用。
本文将简要介绍中心差分的概念和应用，以及动力学方程的概述，并重点探讨 MATLAB 在这两个方面的应用。
二、中心差分的概念和应用中心差分是一种常用的离散差分方法，主要用于离散信号的处理和分析。
中心差分的基本思想是将信号的值按照一定的间隔进行离散，然后用差分代替微分，从而得到离散信号的差分方程。
在实际应用中，中心差分广泛应用于信号处理、图像处理、数据分析等领域。
三、动力学方程的概述动力学方程是描述物体运动规律的数学方程，它包含了物体的质量、速度、加速度等物理量。
动力学方程在物理学、力学、航空航天等领域有着广泛的应用。
根据物体的运动状态和受力情况，动力学方程可以分为牛顿力学方程、拉格朗日方程、哈密顿方程等。
四、MATLAB 在中心差分和动力学方程中的应用MATLAB 是一种功能强大的数学软件，它可以用于解决各种数学问题和工程问题。
在中心差分方面，MATLAB 提供了一系列的差分运算函数，如 diff、gradient 等，用户可以利用这些函数轻松地实现中心差分。
在动力学方程方面，MATLAB 提供了 ODE45、ODE23 等求解器，用户可以利用这些求解器求解各种动力学方程。
五、结论中心差分和动力学方程在科学研究和工程技术中具有重要意义，而MATLAB 作为一种功能强大的数学软件，为解决这两方面的问题提供了方便。
一维中心差分法计算代码1. 引言1.1 一维中心差分法概述一维中心差分法是一种常用的数值计算方法，通常用于求解一维偏微分方程的数值解。
这种方法利用了函数在某一点的导数可以通过函数在该点附近的取值来近似表示的特性。
通过将函数在一个点的导数表示为该点附近两个点的函数值的线性组合，可以得到一维中心差分的计算公式。
一维中心差分法的基本思想是通过离散化空间，将空间区域划分为一系列小区间，然后在每个小区间上使用差分公式来逼近偏微分方程的微分操作。
中心差分是一种常用的差分方式，通过取该点周围两个点的函数值的平均值来表示该点的导数。
在一维中心差分法中，我们可以通过迭代更新的方式，逐步求解出整个空间区域的数值解。
一维中心差分法在数值求解偏微分方程时具有一定的优势，它能够较为准确地近似实际解，并且在计算速度上具有一定的优势。
该方法也存在一些缺点，比如对于非线性问题的数值求解较为困难，且在边界处需要特别处理。
综合考虑其优缺点，一维中心差分法仍然是一个比较常用且有效的数值计算方法。
2. 正文2.1 一维中心差分法计算代码实现步骤一维中心差分法是一种常用的数值计算方法，用于求解一维偏微分方程的数值解。
在本节中，我们将介绍一维中心差分法的计算代码实现步骤，以帮助读者更好地理解这一数值计算方法的具体实现过程。
我们需要定义一个一维数组来存储我们要计算的函数值。
假设我们要求解的一维偏微分方程为u_t = ku_xx，在这里u_t表示u关于时间的导数，k为常数，u_xx表示u关于空间的二阶导数。
我们可以将空间离散化为N个网格点，将时间离散化为M个时间步长，用一个N\*M的数组来存储u在各个网格点上的值。
我们需要初始化u在初始时刻t=0的值。
这可以通过给定的初始条件来完成，比如u(x,0) = f(x)，其中f(x)为已知函数。
然后，我们可以开始计算中心差分来逼近u_xx，即利用有限差分近似u关于空间的二阶导数。
具体地，假设我们在某一时间步t_n和位置x_i处，我们有如下的中心差分公式：∆u_i^n = (u_{i-1}^n - 2u_i^n + u_{i+1}^n) / ∆x^2其中∆u_i^n表示在位置x_i处的u在时间步t_n的二阶空间导数近似值，∆x为空间步长。
中心差分法计算单自由度体系动力反应1，程序说明中心差分法基于有限元查分代替位移对时间的求导（即速度和加速度）。
如果采用等时间步长，i t t ∆=∆，则速度和加速度的中心查分近似为.11..11222i i i i i u u u t u u u u t+-+--=∆-+=∆体系的运动方程为...()()()()mu t cu t ku t P t ++=联立以上三式，得1111222i i i i i i i u u u u u m c ku P t t+-+--+-++=∆∆上式中，假设i u 和1i u -是已知的，即i t 和i t 以前时刻的运动已知，则可以把已知项移到方程的右边，整理得11222222i i i i m c m mc u P k u u t t t tt +-⎛⎫⎛⎫⎛⎫+=---- ⎪ ⎪ ⎪∆∆∆∆∆⎝⎭⎝⎭⎝⎭这样，就可以计算体系任意时刻的位移，速度和加速度。
2，程序框图3，程序清单%计算等效刚度和中心差分法计算公式中的系数clear,m=17.5e3;k=875500;c=35000;aa=input(' 请选择时间步长 1 or 2 or 3 \n 1：dt=0.02 ; 2 ：dt=0.3 ; 3: dt=其它\n');if aa==1dt=0.02;endif aa==2dt=0.3;endif aa==3dt=input('请输入时间步长\n dt= ')endt=0:dt:1.2;n=fix(1.2/dt+1);kr=m/(dt * dt) + c/(2 * dt); a=k-2 * m/(dt * dt); b=m/(dt * dt)-c/(2 * dt);%求力pp1=0:1:40;p2=39:-1:0;one=ones(1,40);p3=(one<0);p=1000*[p1,p2,p3];for i=1:nif t(i)<=0.4,p(i)=100000*t(i);endif t(i)>0.4&&t(i)<=0.8,p(i)=80000-100000*t(i);endif t(i)>0.8,p(i)=0;endend%以下求位移uu(1,1)=0;u(1,2)=0;for i=2:(n-1)u(1,i+1) =( p(1,i) - a * u(1,i) - b * u(1,i-1))/kr;end%以下求速度ui，加速度uiiui(1,1)=0;uii(1,1)=0;for i=2:(n-1)ui(1,i)=(u(1,i+1)-u(1,i-1))/(2 * dt);uii(1,i)=(u(1,i+1)-2 * u(1,i)+u(1,i-1))/(dt*dt);endui(1,n)=2*ui(1,n-1)-ui(1,n-2);uii(1,n)=2*uii(1,n-1)-uii(1,n-2);figure(1)plot(t,u,'b'),xlabel('时间（s）'),ylabel('位移（m）'),title('位移时程曲线 '); figure(2)plot(t,ui,'r'),xlabel('时间（s）'),ylabel('速度（m/s）'),title('速度时程曲线 ');figure(3)plot(t,uii,'k'),xlabel('时间（s ）'),ylabel('加速度(m/s^2)'),title('加速度时程曲线 ');4，输入数据体系质量 m=17.5t刚度 k=875.5kN/m 阻尼 c=35kNs/m 荷载形式{100000(0.4)80000100000(0.40.8)0(0.8)t t s t s t s t s i p <-≤≤>=5，计算结果dt=0.02s 计算结果0.20.40.60.811.21.4-0.04-0.03-0.02-0.0100.010.020.030.040.050.06时间（s ）位移（m ）位移时程曲线0.20.40.60.811.21.4-0.3-0.25-0.2-0.15-0.1-0.0500.050.10.150.2时间（s ）速度（m /s ）00.20.40.60.81 1.2 1.4-2-1.5-1-0.500.511.52时间（s ）加速度(m /s 2)加速度时程曲线dt=0.3s 计算结果00.20.40.60.81 1.2 1.4-0.15-0.1-0.0500.050.10.150.20.25时间（s ）位移（m ）00.20.40.60.81 1.2 1.4-0.3-0.2-0.100.10.20.30.40.50.6时间（s ）速度（m /s ）速度时程曲线00.20.40.60.81 1.2 1.4-55101520时间（s ）加速度(m /s 2)6，稳定性 结构自振周期17.072/0.89n n n rads T sωπω-=====由于0.890.020.283.14nT t s s π∆=<== 所以选t=0.02s 满足稳定性条件。
中心差分法计算程序编程姓名：张泽伟 学号： 电话：一、中心差分法程序原理说明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 的值。
中心差分法计算程序编程具体而言，中心差分法首先选择一个较小的步长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 的值。
第1篇一、实验目的本次实验旨在通过数值分析的方法，对工程实际问题进行建模、求解和分析。
通过学习数值方法的基本原理和算法，提高解决实际工程问题的能力。
二、实验内容1. 线性方程组的求解2. 矩阵特征值与特征向量的计算3. 函数插值与曲线拟合4. 数值微分与积分三、实验步骤1. 线性方程组的求解（1）编写程序实现高斯消元法、克劳斯消元法和列主元素法（2）设计输入界面，用户输入增广矩阵的行和列，填写系数及常数项（3）分别运用三种方法求解线性方程组，比较求解结果的正确性、数值稳定性和计算效率2. 矩阵特征值与特征向量的计算（1）编写程序实现幂法、QR算法和逆幂法（2）设计输入界面，用户输入矩阵的行和列，填写矩阵元素（3）分别运用三种方法计算矩阵的特征值与特征向量，比较求解结果的准确性和计算效率3. 函数插值与曲线拟合（1）编写程序实现拉格朗日插值、牛顿插值和样条插值（2）设计输入界面，用户输入函数的自变量和函数值，选择插值方法（3）分别运用三种方法进行函数插值，比较插值结果的准确性和光滑性4. 数值微分与积分（1）编写程序实现有限差分法、龙格-库塔法和辛普森法（2）设计输入界面，用户输入函数的导数或积分的上下限，选择数值方法（3）分别运用三种方法进行数值微分和积分，比较求解结果的准确性和计算效率四、实验结果与分析1. 线性方程组的求解通过实验，我们发现列主元素法在求解线性方程组时具有较好的数值稳定性，计算效率也较高。
而高斯消元法和克劳斯消元法在处理大型稀疏矩阵时存在一定的困难。
2. 矩阵特征值与特征向量的计算实验结果表明，QR算法和逆幂法在计算矩阵特征值与特征向量时具有较高的准确性和计算效率。
幂法在处理大型稀疏矩阵时表现出较好的性能。
3. 函数插值与曲线拟合在函数插值和曲线拟合实验中，样条插值方法具有较好的准确性和光滑性。
拉格朗日插值和牛顿插值方法在处理简单函数时表现良好，但在处理复杂函数时可能存在精度问题。
中心差分格式1、考虑问题考虑二阶常微分方程边值问题：f qu dxu d Lu =+-=22 （1） βα==)(,)(b u a u其中q ，f 为[a,b]上的连续函数，βα,为常数。
2、网格剖分与差分格式将区间[a,b]分成N 等分,分点为N i ih a x i ,,1,0,⋅⋅⋅=+=,h=(b-a)/N,于是我们得到区间I=[a,b]的网格剖分，i x 为网格节点，h 为步长。
差分格式为：.,,1,,2,120211βα==-⋅⋅⋅==++--=-+N i i i i i i i h u u N i f u q hu u u u L3、截断误差将方程（1）在节点离散化，由泰勒公式展开得)()(12)()()(2)(344222211h dx x u d h dx x u d h x u x u x u ii i i i O +⎥⎦⎤⎢⎣⎡+⎥⎦⎤⎢⎣⎡=+--+ 所以截断误差为)()(12)(3442h dx x u d h u R ii O +⎥⎦⎤⎢⎣⎡-= 4、数值例子x x q e x u x sin 1)()(+== 其中[]1,0∈x5、求解 由f qu dxu d Lu =+-=22，且已知 x x q e x u xsin 1)()(+== 可得x e x f x sin )(=将向量式的差分格式用矩阵形式表示出来，得到矩阵形式为⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎣⎡+--+--+-2122212112112h q h q h q N⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡-121N u u u =⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎣⎡++-βα122212N f h f h f h 系数矩阵A=⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎣⎡+--+--+-2122212112112h q h q h q N ，我们求出矩阵A 极其逆便可求得u （x ）的数值解。
6、参考文献《偏微分方程数值解法》李荣华 高等教育出版社《科学计算中的有限差分法》《MATLAB 程序设计教程》刘卫国 中国水利水电出版社附件1、程序流程图开始计算各初始条件构造系数矩阵d计算矩阵d的逆D构造格式右端矩阵B计算数值解u=D*B输出u的数值解结束2 、程序代码%zah.m%u(x)=exp(x), q(x)=1+sinx, f(x)=exp(x)*sinx 0<x<1a=input('input a:');%输入区间条件ab=input('input b:');%输入区间条件bN=input('input N:');%输入等分数Nh=(b-a)/N;%步长for i=1:N-1x(i)=a+i*h;%网格节点q(i)=1+sin(x(i));%连续函数q网格节点上各值f(i)=exp(x(i))*sin(x(i));%连续函数f网格节点上各值endc1=linspace(2+q(1)*h*h,2+q(N-1)*h*h,N-1); %系数矩阵主对角线上各值c2=linspace(-1,-1,N-2); %系数矩阵次对角线上各值d1=diag(c1,0);d2=diag(c2,-1);d3=diag(c2,1);d=d1+d2+d3 %构造系数矩阵D=inv(d); %系数矩阵的逆b1=linspace(h*h*f(1),h*h*f(N-1),N-1);for i=1:N-1if i==1b2(i)=exp(a);elseif i==N-1b2(i)=exp(b);elseb2(i)=0;endendendB1=b1+b2;%构造差分格式等式右边的向量B=B1' %差分格式等式右边向量的转置u=D*B %输出u（x）数值解输出结果：>> clear>> zahinput a:0input b:1input N:10d =2.0110 -1.0000 0 0 0 0 0 0 0 -1.0000 2.0119 -1.0000 0 0 0 0 0 0 0 -1.0000 2.0127 -1.0000 0 0 0 0 0 0 0 -1.0000 2.0136 -1.0000 0 0 0 0 0 0 0 -1.0000 2.0144 -1.0000 0 0 0 0 0 0 0 -1.0000 2.0153 -1.0000 0 0 0 0 0 0 0 -1.0000 2.0161 -1.0000 0 0 0 0 0 0 0 -1.0000 2.0170 -1.0000 0 0 0 0 0 0 0 -1.0000 2.0178B =1.00110.00340.00560.00790.01020.01250.01470.01702.7375u =1.11171.23451.36851.51431.67271.84512.03312.23932.4664。
中心差分法的基本理论与程序设计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)按照相同的步骤，所得结果如下：再计算下去，位移将继续增大，这是不稳定的典型表现。
其原因是在条件稳定的可能得到有意义的结果。
6.2.2振型叠加法（精确解）采用振型叠加法对上述问题进行计算，可以得到该问题的精确解。
首先应求解的广义特征值问题为：按照一般的线性代数方法可以得到上式的解答为：利用上式，3个互不耦合的运动方程，即：3+()i t=利用无阻尼情形的杜哈美积分公式，可以得到上述方程的解答为：最后利用振型叠加得到系统的位移为：根据上式计算得到的每一时间步长的位移值是系统响应的精确解，如表2 和表3所示。
(1)表20.00 0.00 0.01 0.04 0.14 0.37 0.79 1.45 2.32 3.34 0.00 0.04 0.21 0.59 1.25 2.21 3.38 4.63 5.80 6.76 0.391.452.924.485.816.747.297.607.928.46(2)表33.89 3.46 -1.19 2.25 1.83 -2.40 1.78 2.25 -1.23 3.40 6.75 6.76 0.00 6.73 6.77 0.00 6.72 6.79 0.00 6.70 8.178.410.608.979.241.209.199.050.618.366.3 程序计算(1)程序的计算过程如图2所示，具体输出结果如表4所示。
表4 3100.363t T ∆==数值解时间t ∆ 2t ∆ 3t ∆ 4t ∆ 5t ∆ 6t ∆ 7t ∆ 8t ∆ 9t ∆ 10t ∆ 1a0.00 0.00 0.01 0.03 0.13 0.36 0.79 1.46 2.37 3.42 2a 0.00 0.04 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 3100.363t T ∆==时的程序计算过程及计算结果(2) 时间步长3518.14t T ∆==当3518.14t T ∆==时，计算得到的起步条件为[]00987.179Tt a -∆=该程序的计算过程如图3所示，具体输出结果如表5所示。
表5 3518.14t T ∆==数值解时间t ∆ 2t ∆ 3t ∆4t ∆ 5t ∆1a0.00 0.0077.1310⨯111.2410-⨯ 141.5910⨯2a 0.0052.1710⨯ 82.3610-⨯ 112.3510⨯ 142.3210-⨯3a29.8710⨯56.4610-⨯85.6610⨯115.2710-⨯ 145.0110⨯图3 3518.14t T ∆==时的程序计算过程及计算结果6.4 结果讨论(1) 时间步长3100.363t T ∆==对于3100.363t T ∆==的情况，三者的比较如表6和图4所示。
表6 3100.363t T ∆==理论解，精确解和数值解比较结果理论解精确解数值解t(s) a1(m) a2(m) a3(m) a1(m) a2(m) a3(m) a1(m) a2(m) a3(m)0.726 0 0 0.4 0 0 0.39 0 0 0.41.089 0 0.03 1.48 0 0.04 1.45 0 0.04 1.48 1.452 0 0.192.97 0.01 0.21 2.92 0.01 0.19 2.971.815 0.03 0.58 4.52 0.04 0.59 4.48 0.03 0.58 4.522.178 0.13 1.26 5.82 0.14 1.25 5.81 0.13 1.26 5.82 2.541 0.36 2.24 6.71 0.37 2.21 6.74 0.36 2.24 6.712.904 0.793.43 7.22 0.79 3.38 7.29 0.79 3.43 7.223.267 1.464.69 7.51 1.45 4.63 7.6 1.46 4.69 7.51 3.63 2.375.84 7.85 2.32 5.8 7.92 2.37 5.84 7.85 3.993 3.426.77 8.45 3.34 6.76 8.46 3.42 6.77 8.45(a)(b)图4，不论用理论计算还是程序计算，都与精确解吻合的非常好。
说明了该程序用于计算实际问题的准确性和有效性，能够满足精度要求。
(2)程序计算，所得到的位移都会无限增大，说明此时的算法是不稳定的，与精确解并没有可比性。
这也可以证明中心差分法是条件稳定算法，在利用其求解具体问否则算法将是不稳定的。
为了保证解的稳定性，中心差分法的步长必须满足下列条件：估计得到。
7总结本文中的程序采用了C语言（部分C++语言）实现了用于求解动力学问题的中心差分法，可以求解得到运动方程的解答，包括位移，速度和加速度。
计算简便且在算法稳定的条件下，精度较高。
中心差分法是对运动方程不进行方程形式的变换而直接进行逐步数值积分，属于直接积分法。
中心差分法的基本思路是用有限差分代替位移对时间的求导，将运动方程中的速度和加速度用位移的某种组合表示，然后将常微分方程组的求解问题转换为代数方程组的求解问题，并假设在每个小的时间间隔内满足运动方程，则可以求得每个时间间隔的递推公式，进而求得整个时程的反应。
中心差分法的基本算法步骤为：输入刚度矩阵，质量矩阵和阻尼矩阵，并且给定初始时刻的位移，速度，加速度（可通过计算得到）；计算积分常数；计算起步条件；形成有效质量矩阵，并对其进行三角分解；对于求解。
本文通过求解一个三自由度系统的运动方程验证了该程序的正确性。
首先通过中心差分法，振型叠加法和该程序求解得到了问题的理论解，精确解和数值解。
解和数值解吻合的非常好，证明该算法的合理性与准确性；理论解和数值解得到的位移结果都会无限增大，说明此时的算法是不稳定的，与精确解并没有可比性，这也间接证明中心差分法是条件稳定算法。
8心得体会我认为在编写程序的过程中，设计算法是核心部分，需要对基础理论有很深刻的理解，清楚求解方法的每一个步骤，同时还要考虑程序语言的实现。
在得到算法步骤以后，需要将其翻译成计算机程序设计语言，然后进行编辑，调试和运行，得到运行结果。
得到运行结果并不意味着程序正确，还需要算例对其进行验证分析。
附录// 中心差分法程序#include <stdio.h>#include <iostream.h>#include <iomanip.h>#include <fstream.h>using namespace std;//调用子程序void ArrayLU(double**, double**, double**, int);void Inverse(double**, double**, int);void ArrayMVector(double**, double*, double*, int);void LUSolve(double**, double**, double*, double*,int);void Centre(double**, double**, double**, double*, double**, double**, double**, int, int, double);int main(){//设置矩阵维度，时间步，时间步长int num=3;int timeStep=11;double dt=0.363; // double dt=18.14;//开辟矩阵空间double **M = new double*[num];for(int i=0; i<num; i++)M[i] = new double[num];double **C = new double*[num];for(int i=0; i<num; i++)C[i] = new double[num];double **K = new double*[num];for(int i=0; i<num; i++)K[i] = new double[num];double **X = new double*[timeStep+1];for(int i=0; i<timeStep+1; i++)X[i] = new double[num];double **V = new double*[timeStep+1];for(int i=0; i<timeStep+1; i++)V[i] = new double[num];double **A = new double*[timeStep+1];for(int i=0; i<timeStep+1; i++)A[i] = new double[num];double *P = new double[num];//设置计算参数与初始条件X[1][0] = 0; X[1][1] = 0; X[1][2] = 0;V[1][0] = 0; V[1][1] = 0; V[1][2] = 0;A[1][0] = 0; A[1][1] = 0; A[1][2] = 6;M[0][0] = 1; M[0][1] = 0; M[0][2] = 0;M[1][0] = 0; M[1][1] = 3; M[1][2] = 0;M[2][0] = 0; M[2][1] = 0; M[2][2] = 1;C[0][0] = 0; C[0][1] = 0; C[0][2] = 0;C[1][0] = 0; C[1][1] = 0; C[1][2] = 0;C[2][0] = 0; C[2][1] = 0; C[2][2] = 0;K[0][0] = 2; K[0][1] = -1; K[0][2] = 0;K[1][0] = -1; K[1][1] = 4; K[1][2] = -2;K[2][0] = 0; K[2][1] = -2; K[2][2] = 2;P[0]=0; P[1]=0; P[2]=6;//开始计算，调用子程序CentreCentre(M,C,K,P,X,V,A,num,timeStep,dt);//输出计算结果ofstream out("out.txt"); //输出结果文件out.txtif(!out.is_open())return;for(int i=0; i<=timeStep; i++){cout<<i*dt<<" s:"<<endl; //屏幕显示out<< i*dt<<" s:"<<endl; //将结果写入结果文件for(int j=0; j<num; j++) {cout << setw(20) << X[i][j] << setw(20) << V[i][j] << setw(20) << A[i][j] << endl; out<<setw(20) << X[i][j] << setw(20) << V[i][j] << setw(20) << A[i][j] << endl;} }out.close();return 0;}//子程序Centrevoid Centre(double**M, double**C, double**K, double *P, double**X, double**V, double**A, int size, int timeStep, double dt){double c0,c1,c2,c3;double *PT=new double[size];double *temp=new double[size];double *tran=new double[size];double **MM=new double*[size];for(int i=0; i<size; i++)MM[i]=new double[size];double **L = new double* [size];for(int i=0; i<size; i++)L[i] = new double [size];double **U = new double* [size];for(int i=0; i<size; i++)U[i] = new double [size];//计算积分常数c0=1/(dt*dt);c1=1/(2*dt);c2=2/(dt*dt);c3=(dt*dt)/2;//X[0][0]=X[1][0]-dt*V[1][0]+c3*A[1][0];X[0][1]=X[1][1]-dt*V[1][1]+c3*A[1][1];X[0][2]=X[1][2]-dt*V[1][2]+c3*A[1][2];//形成有效质量矩阵for(int i=0; i<size; i++){for(int j=0; j<size; j++)MM[i][j]=c0*M[i][j]+c1*C[i][j];}//调用子程序ArrayLU对有效质量矩阵进行LU分解ArrayLU(MM,L,U,size);//计算有效载荷for(int ti=1; ti<=timeStep-1; ti++){for(int i=0; i<size; i++)temp[i] = c2*X[ti][i]-c0*X[ti-1][i];//调用ArrayMVector，得到矩阵和向量乘积ArrayMVector(M,temp,tran,size);for(int i=0; i<size; i++)PT[i] = P[i]+tran[i];for(int i=0; i<size; i++)temp[i] = c1*X[ti-1][i];ArrayMVector(C,temp,tran,size);for(int i=0; i<size; i++)PT[i] = PT[i]+tran[i];for(int i=0; i<size; i++)temp[i] = X[ti][i];ArrayMVector(K,temp,tran,size);for(int i=0; i<size; i++)PT[i] = PT[i]-tran[i];//LUSolve(L,U,PT,temp,size);//输出结果：位移，速度，加速度for(int i=0; i<size; i++){X[ti+1][i] = temp[i];A[ti][i] = c0*(X[ti-1][i]-2*X[ti][i]+X[ti+1][i]);V[ti][i] = c1*(-X[ti-1][i]+X[ti+1][i]);}}//删除不用数组for(int i=0; i<size; i++)delete [] MM[i];delete []MM;for(int i=0; i<size; i++)delete [] L[i];delete []L;for(int i=0; i<size; i++)delete [] U[i];delete []U;delete []PT;delete []temp;delete []tran;}//子程序ArrayLU，对有效质量矩阵进行LU分解void ArrayLU(double**A, double**L, double**U, int size) {double scale;double **B=new double*[size];for(int i=0; i<size; i++)B[i]=new double [size];for(int i=0; i<size; i++){for(int j=0; j<size; j++){B[i][j] = 0;L[i][j] = 0;U[i][j] = 0;}B[i][i]=1;}for(int i=0; i<size-1; i++){if(A[i][i]==0)cout<<"the Array is singular"<<endl;for(int j=i+1; j<size; j++){scale = A[j][i]/A[i][i];A[j][k] = A[j][k]-A[i][k]*scale;for(int k=0; k<size; k++)B[j][k] = B[j][k]-B[i][k]*scale;}}//调用子程序Inverse，对矩阵进行转置Inverse(B, L, size);for(int i=0; i<size; i++){for(int j=i; j<size; j++)U[i][j] = A[i][j];}for(int i=0; i<size; i++)delete[] B[i];delete[] B;}//子程序Inverse，对矩阵进行转置void Inverse(double**A, double**B, int size) {double scale;for(int i=0; i<size; i++)B[i][i] = 1;for(int i=0; i<size-1; i++){if(A[i][i]==0)cout << "the Array is singular" << endl;scale = A[j][i]/A[i][i];for(int k=i; k<size; k++)A[j][k] = A[j][k]-A[i][k]*scale;for(int k=0; k<size; k++)B[j][k] = B[j][k]-B[i][k]*scale;}}for(int i=size-1; i>0; i--){for(int k=0; k<size; k++)B[i][k] = B[i][k]/A[i][i];A[i][i] = 1;for(int j=i-1; j>=0; j--){scale = A[j][i]/A[i][i];A[j][i] = 0;for(int k=0; k<size; k++)B[j][k] = B[j][k]-B[i][k]*scale;}}for(int k=0; k<size; k++)B[0][k] = B[0][k]/A[0][0];}//子程序ArrayMVector，矩阵和向量的乘法void ArrayMVector(double **A, double *X, double *B, int size) {for(int i=0; i<size; i++){B[i] = 0;for(int j=0; j<size; j++)B[i] = B[i] + A[i][j]*X[j];}}//子程序LUSolvevoid LUSolve(double**L, double**U, double*P, double*X, int size) {double sum;double* Y = new double[size];Y[0]=P[0]/L[0][0];for(int i=1; i<size; i++){sum = 0;for(int j=0; j<i; j++)sum+=L[i][j]*Y[j];Y[i]=(P[i]-sum)/L[i][i];}X[size-1] = Y[size-1]/U[size-1][size-1];for(int i=size-2; i>=0; i--){sum=0;for(int j=i+1; j<size; j++)sum = sum + U[i][j]*X[j];X[i]=(Y[i]-sum)/U[i][i];}delete []Y; }。
