中心差分法的基本理论与程序设计

中心差分法的基本理论与程序设计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; }。