中心差分法計(jì)算程序編程

1、合用標(biāo)準(zhǔn)文案中心差分法計(jì)算程序編程姓名：張澤偉學(xué)號(hào)：電話：一、中心差分法程序原理說明1.1中心差分法思路中心差分法的基本思路：是將運(yùn)方程中的速度向量和加速度向量用位移的某種合來表示，將微分方程的求解化代數(shù)方程的求解，并在區(qū)內(nèi)求得每個(gè)渺小區(qū)的推公式，而求得整個(gè)程的反。1.2中心差分法原理中心差分法只在相隔t一些失散的區(qū)內(nèi)足運(yùn)方程，其基于有限差分代替位移的求（即速度和加速度），若是采用等步，tit，速度與加速度的中心差分近似：uui1ui12t(a)uui12uiui1t2(b)而失散點(diǎn)的運(yùn)uiu(ti),uiu(ti),uiu(ti)(i0,1,2,3,)由系統(tǒng)運(yùn)方程：mu(ti)cu(t)ku

2、(t)0(c)將速度和加速度的差分近似公式（a）和式（b）代入式（c）可以獲取ti刻的運(yùn)方程：mui12uiui1cui1ui1kui0t22t（d）在（d）式中，假ui和ui1是已知的，即在ti及ti以前刻的運(yùn)已知，可以把已知移到方程的右，整理獲?。?mc)ui1(k2m)ui(mc)ui1(e)t22tt2t22t優(yōu)秀文檔合用標(biāo)準(zhǔn)文案由式（e）就可以依照ti及ti以前時(shí)刻的運(yùn)動(dòng)，求得ti1時(shí)刻的運(yùn)動(dòng)，若是需要可以用式（a）和式（b）求得系統(tǒng)的速度和加速度。1.3初始條件轉(zhuǎn)變假設(shè)給定的初始條件為u0u(0),u0u(0),（g）由式（g）確定u1。在零時(shí)刻速度和加速度的中心差分公式為：u1u

3、1u02t（h）u12u0u1u0t2（i）u1u0tu0t2u0（j）將式（i）消去u1得：2而零時(shí)刻的加速度值u0可以用t0時(shí)的運(yùn)動(dòng)方程mu0cu0ku00確定u01(cu0ku0)即m（k）這樣就可以依照初始條件u0,u0和初始荷載P0，就可以依照上式確定u1的值。1.4中心差分法編程思路基本數(shù)據(jù)準(zhǔn)備和初始條件計(jì)算：u01(cu0ku0)mu1u0tu0t2u02計(jì)算等效剛度和中心差分計(jì)算公式中的相關(guān)系數(shù)：mck22tt優(yōu)秀文檔合用標(biāo)準(zhǔn)文案ak2mt2mcbt22t依照ti及ti以前時(shí)刻的運(yùn)動(dòng)，計(jì)算ti1時(shí)刻的運(yùn)動(dòng)：Pauibui1ui1Pkui1ui1u2tui12uiui1ut2下一

4、步計(jì)算用i+1代替i，對(duì)于線彈性結(jié)構(gòu)系統(tǒng)，重復(fù)第3步，對(duì)于非線性結(jié)構(gòu)系統(tǒng)，重復(fù)第2步和第3步。1.5中心差分法牢固條件以上為中心差分法漸漸計(jì)算公式，其擁有2階精度，即誤差0(t2)；而且為t有條件牢固，牢固條件為：二、程序框圖Tn依照中心差分法的原理，可以得出本程序的主要程序思想，以下面框圖的形式顯現(xiàn)出來：編寫回調(diào)函數(shù)調(diào)用輸入數(shù)據(jù)計(jì)算初始加速度和位移前一時(shí)刻值開始循環(huán)體（步數(shù)控制）計(jì)算等效剛度和系數(shù)a、b計(jì)算位移、速度、加速度并保存優(yōu)秀文檔合用標(biāo)準(zhǔn)文案解析牢固條件繪圖三、程序清單%m,k,c分別為質(zhì)量、剛度、阻尼%p0,dt,t分別為外荷載幅值、時(shí)間步距、總時(shí)間%u0,v0為初始條件初位移和初

5、速度%u,v,ac分別為位移、速度、加速度反應(yīng)ek=等效剛度；p=荷載；ep=等效荷載%定義矩陣X0=input(請(qǐng)按格式溫序次輸入初始矩陣，如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);%分別取出其中的參數(shù)：v0=X0(1,5);t=X0(1,6);P0=X0(1,7);dt=X0(1,8)t=0:dt:t;%將時(shí)間分步，采用等時(shí)間步長；mm,nn=size(t);%計(jì)算t的向量長度，得出步數(shù)；u=zeros(size(t);%設(shè)定儲(chǔ)藏u的矩陣；v=zeros(size(t);%設(shè)定儲(chǔ)藏v的矩陣；ac=zer

6、os(size(t);%設(shè)定儲(chǔ)藏ac的矩陣；u(:,2)=u0;%賦值向量第2項(xiàng)為u0；v(:,2)=v0;%賦值向量第2項(xiàng)為v0；ac(:,2)=(P0-c*v(:,2)-k*u(:,2)/m;%求出初始加速度ac0；u(:,1)=u(:,2)-dt*v(:,2)+(dt)2)*ac(:,2)/2;%計(jì)算初始條件u-1項(xiàng)；ek=m/(dt2)+c/(2*dt);%計(jì)算等效剛度；a=k-(2*m)/(dt2);b=m/(dt2)-c/(2*dt);%計(jì)算方程系數(shù)；p(:,2)=P0*sin(0);%給出初始荷載條件；ep(:,2)=p(:,2)-a*u(:,2)-b*u(:,1);%計(jì)算初始等

7、效荷載；u(:,3)=ep(:,2)/ek;%計(jì)算位移u1=u(:,3)fori=3:nn%從第二項(xiàng)開始進(jìn)行中心差分法計(jì)算；p(:,i)=P0*sin(.5*pi*(i-2)*dt);%給出荷載條件，依照簡諧荷載計(jì)算；ep(:,i)=p(:,i)-a*u(:,i)-b*u(:,i-1);%計(jì)算等效荷載；%-得出所需要結(jié)果-%u(:,i+1)=ep(:,i)/ek;%計(jì)算位移量；v(:,i)=(u(:,i+1)-u(:,i-1)/(2*dt);%計(jì)算速胸襟；ac(:,i)=(u(:,i+1)-2*u(:,i)+u(:,i-1)/(dt2);%計(jì)算加速胸襟；endt=t(:,1:end-1);u=

8、u(:,2:end-1);v=v(:,2:end);ac=ac(:,2:end);p=p(:,2:end);ep=ep(:,2:end);%-繪制位移、速度、加速度時(shí)程曲線-%plot(t,u,b-o),holdon,plot(t,v,g-p),holdon,plot(t,ac,r:x),gridon,xlabel(時(shí)間(s),ylabel(位移(m)速度(m/s)加速度(m/s2),title(頂層u,v,ac的時(shí)程曲線);subplot(3,1,1),plot(t,u,b-),grid,xlabel(時(shí)間(s),ylabel(位移(m),title(位移u的時(shí)程曲線);legend(位移u

9、)優(yōu)秀文檔合用標(biāo)準(zhǔn)文案subplot(3,1,2),plot(t,v,k),grid,xlabel(時(shí)間(s),ylabel(速度(m/s),title(速度v的時(shí)程曲線);legend(速度v)subplot(3,1,3),plot(t,ac,r),grid,xlabel(時(shí)間(s),ylabel(加速度(m/s2),title(加速度ac的時(shí)程曲線);legend(加速度ac)四、輸入數(shù)據(jù)本程序采用單自由度系統(tǒng)進(jìn)行計(jì)算，主要已知參數(shù)信息以下：其質(zhì)量M=9240kg、剛度K1460KN/m、阻尼系數(shù)C6.41kNs/m，對(duì)結(jié)構(gòu)施加動(dòng)力荷載P73000sin(0.5t)N，結(jié)構(gòu)周期T=0.05s，初始位移u00.05m，初始速度v00m/s，假設(shè)結(jié)構(gòu)處于線彈性狀態(tài)。由中心差分法可知，要使計(jì)算結(jié)果牢固且不發(fā)散，需滿足：時(shí)間步長tTn0.159s，本例分別取時(shí)間步長為0.1s、0.15s、優(yōu)秀文檔合用標(biāo)準(zhǔn)文案當(dāng)dt=0.15s時(shí)：當(dāng)dt=0.17s時(shí)：優(yōu)秀文檔合用標(biāo)準(zhǔn)文案當(dāng)dt=0.2s時(shí)：六、結(jié)果牢固性解析由以上時(shí)程圖可以獲適當(dāng)t0.1，0.15時(shí)漸漸計(jì)算結(jié)果給出的結(jié)構(gòu)運(yùn)動(dòng)趨向收斂的，即計(jì)