中心差分法计算程序编程

《中心差分法计算程序编程》由会员上传分享，免费在线阅读，更多相关内容在工程资料-天天文库。

1、实用标准文案中心差分法计算程序编程姓名：张泽伟学号：电话：一、中心差分法程序原理说明1.1中心差分法思路中心差分法的基本思路：是将运动方程中的速度向量和加速度向量用位移的某种组合来表示，将微分方程组的求解问题转化为代数方程组的求解问题，并在时间区间内求得每个微小时间区间的递推公式，进而求得整个时程的反应。1.2中心差分法原理中心差分法只在相隔一些离散的时间区间内满足运动方程，其基于有限差分代替位移对时间的求导（即速度和加速度），如果采用等时间步长，，则速度与加速度的中心差分近似为：(a)(b)而离散时间点的运动为(0,1,2,3,……)由体系运动方

2、程为：(c)将速度和加速度的差分近似公式（a）和式（b）代入式（c）可以得到时刻的运动方程：（d）在（d）式中，假设和是已知的，即在及以前时刻的运动已知，则可以把已知项移到方程的右边，整理得到：(e)精彩文档实用标准文案由式（e）就可以根据及以前时刻的运动，求得时刻的运动，如果需要可以用式（a）和式（b）求得体系的速度和加速度。1.3初始条件转化假设给定的初始条件为（g）由式（g）确定。在零时刻速度和加速度的中心差分公式为：（h）`（i）将式（i）消去得：（j）而零时刻的加速度值可以用t＝0时的运动方程确定即（k）这样就可以根据初始条件和初始荷载，

3、就可以根据上式确定的值。1.4中心差分法编程思路①基本数据准备和初始条件计算：②计算等效刚度和中心差分计算公式中的相关系数：精彩文档实用标准文案③根据及以前时刻的运动，计算时刻的运动：④下一步计算用i+1代替i，对于线弹性结构体系，重复第3步，对于非线性结构体系，重复第2步和第3步。1.5中心差分法稳定条件以上为中心差分法逐步计算公式，其具有2阶精度，即误差；并且为有条件稳定，稳定条件为：二、程序框图根据中心差分法的原理，可以得出本程序的主要程序思想，以下面框图的形式展示出来：编辑回调函数调用输入数据计算初始加速度和位移前一时刻值开始循环体（步数控

4、制）计算等效刚度和系数a、b计算位移、速度、加速度并保存分析稳定条件精彩文档实用标准文案绘图二、程序清单%m,k,c分别为质量、刚度、阻尼%p0,dt,t分别为外荷载幅值、时间步距、总时间%u0,v0为初始条件初位移和初速度%u,v,ac分别为位移、速度、加速度反应ek=等效刚度；p=荷载；ep=等效荷载%定义矩阵X0=input('请按格式和顺序输入初始矩阵，如X0=[m,k,c,u0,v0,t,P0,dt],m=X0(1,1);k=X0(1,2);c=X0(1,3);u0=X0(1,4);%分别取出其中的参数：v0=X0(1,5);t=X0(1

5、,6);P0=X0(1,7);dt=X0(1,8)t=[0:dt:t];%将时间分步，采用等时间步长；[mm,nn]=size(t);%计算t的向量长度，得出步数；u=zeros(size(t));%设定存储u的矩阵；v=zeros(size(t));%设定存储v的矩阵；ac=zeros(size(t));%设定存储ac的矩阵；u(:,2)=u0;%赋值向量第2项为u0；v(:,2)=v0;%赋值向量第2项为v0；ac(:,2)=(P0-c*v(:,2)-k*u(:,2))/m;%求出初始加速度ac0；u(:,1)=u(:,2)-dt*v(:,2)

6、+((dt)^2)*ac(:,2)/2;%计算初始条件u-1项；ek=m/(dt^2)+c/(2*dt);%计算等效刚度；a=k-(2*m)/(dt^2);b=m/(dt^2)-c/(2*dt);%计算方程系数；p(:,2)=P0*sin(0);%给出初始荷载条件；ep(:,2)=p(:,2)-a*u(:,2)-b*u(:,1);%计算初始等效荷载；u(:,3)=ep(:,2)/ek;%计算位移u1=u(:,3)fori=3:nn%从第二项开始进行中心差分法计算；p(:,i)=P0*sin(.5*pi*(i-2)*dt);%给出荷载条件，按照简谐荷

7、载计算；ep(:,i)=p(:,i)-a*u(:,i)-b*u(:,i-1);%计算等效荷载；%-----------------------得出所需要结果----------------------------------%u(:,i+1)=ep(:,i)/ek;%计算位移量；v(:,i)=(u(:,i+1)-u(:,i-1))/(2*dt);%计算速度量；ac(:,i)=(u(:,i+1)-2*u(:,i)+u(:,i-1))/(dt^2);%计算加速度量；endt=t(:,1:end-1);u=u(:,2:end-1);v=v(:,2:end

8、);ac=ac(:,2:end);p=p(:,2:end);ep=ep(:,2:end);%---------------