M=[2,0;0 1 ];? ?? ?? ?? ?? ?? ?? ? %質量矩陣

K=[6 -2;-2 4];? ?? ?? ?? ?? ?? ? %剛度矩陣

a=0;b=0;

C=a*K+b*M;

dt=0.28;

t=0:dt:2.8;

ft0=zeros(length(K),length(t));

for i=1:length(t)

ft0(1,i)=10;??%在節點4的豎直方向加大小為200N的階躍力

end

dsp=zeros(length(K),length(t));? ?? ?? ?? ?? ?? ?? ?? ?? ?? ?? ?? ?? ???% 位移

vel=zeros(length(K),length(t));? ?? ?? ?? ?? ?? ?? ?? ?? ?? ?? ?? ?? ?? ?? ?% 速度

acc=zeros(length(K),length(t));? ?? ?? ?? ?? ?? ?? ?? ?? ?? ?? ?? ?? ?? ?% 加速度

%--------------------------------------------------------------------------

%??(2) Newmark

alpha=0.3; beta=0.6;? ?? ?? ?? ?? ?? ?? ?? ?? ?? ?? ?? ?? ? % 穩定條件

acc(:,1)=inv(M)*(ft0(:,1)-K*dsp(:,1)-C*vel(:,1));? ?% 計算初始加速度 (t=0)

ekk=K+M/(alpha*dt^2)+C*beta/(alpha*dt);??% 計算有效剛度矩陣

for it=2:length(t)-1? ?? ?? ?? ?? ?? ?? ?? ?? ?? ?? ?? ?? ?? ?? ? % 時間步循環

cfm=dsp(:,it)/(alpha*dt^2)+vel(:,it)/(alpha*dt)+acc(:,it)*(0.5/alpha-1);

cfc=dsp(:,it)*beta/(alpha*dt)+vel(:,it)*(beta/alpha-1)...

+acc(:,it)*(0.5*beta/alpha-1)*dt;

efd=ft0(:,it)+M*cfm+C*cfc;? ?? ?? ?? ?? ?? ?%??計算有效力矢量

dsp(:,it+1)=inv(ekk)*efd;? ?? ?? ?? ?? ?? ?? ?? ?%??t+dt時刻的位移

acc(:,it+1)=(dsp(:,it+1)-dsp(:,it))/(alpha*dt^2)-vel(:,it)/(alpha*dt)...

-acc(:,it)*(0.5/alpha-1);? ?? ?? ?? ???%??t+dt時刻的加速度

vel(:,it+1)=vel(:,it)+acc(:,it)*(1-beta)*dt+acc(:,it+1)*beta*dt; %??t+dt時刻的速度

end

plot(t, dsp(2,:),'-b*')

hold on

xlabel('Time(seconds)')

ylabel(' Vertical displ. (m)')

%??(2) Wilson

%--------------------------------------------------------------------------

theta=1.4;? ?? ?? ?? ?? ?? ?? ?? ?? ?? ?? ?? ?? ?? ?% 穩定條件參數

acc(:,1)=inv(M)*(ft0(:,1)-K*dsp(:,1)-C*vel(:,1)); % 計算初始加速度 (t=0)

ekk=K+M*(6/(theta*dt)^2)+C*(3/(theta*dt));??% 計算有效剛度矩陣

for it=2:length(t)-1

cfm=dsp(:,it)*(6/(theta*dt)^2)+vel(:,it)*(6/(theta*dt))+2*acc(:,it);

cfc=dsp(:,it)*(3/(theta*dt))+2*vel(:,it)+acc(:,it)*(theta*dt/2);

efd=ft0(:,it)+theta*(ft0(:,it+1)-ft0(:,it))+M*cfm+C*cfc; %計算有效力矢量

dtheta=inv(ekk)*efd;? ?? ?? ?? ?? ?? ?? ?? ? %??t+ dt時刻的位移

acc(:,it+1)=(dtheta-dsp(:,it))*(6/(theta^3*dt^2))...

-vel(:,it)*(6/(theta^2*dt))+acc(:,it)*(1-3/theta);??%??t+dt時刻的加速度

vel(:,it+1)=vel(:,it)+acc(:,it+1)*dt/2+acc(:,it)*dt/2; %??t+dt時刻的速度

dsp(:,it+1)=dsp(:,it)+vel(:,it)*dt...

+(acc(:,it+1)+2*acc(:,it))*(dt^2/6);? ?? ?? ?%??t+dt時刻的位移

end

plot(t, dsp(2,:),'-g*')

hold on

xlabel('Time(seconds)')

ylabel(' Vertical displ. (m)')

%--------------------------------------------------------------------------

%??(2) 振型疊加法

%施加不平衡激振力

t=t';

[V,D]=eig(K,M);? ?? ?? ?? ?? ?? ?? ???% 計算特征值和特征向量

[lambda,ki]=sort(diag(D));? ?? ?? ?? ?? ?? ?? ?% 給特征值和特征向量排序

omega=sqrt(lambda);? ?? ?? ?? ?? ?? ?? ?? ?? ?? ?? ?? ?? ?% 角頻率w.

omega1=sqrt(lambda)/(2*pi);? ?? ?? ?? ?? ?? ?? ?? ?? ?% 固有頻率 Hz.

V1=V(:,ki);

h=V1'*M*V1;

Factor=diag(V1'*M*V1);? ?? ?? ?? ?? ?? ?? ?? ?? ???%模態質量

Vnorm=V1*inv(sqrt(diag(Factor)));? ?? ?? ?? ?? ?? ?% 正則化振型向量

VnormK=diag(Vnorm'*K*Vnorm);? ?? ?? ?? ?? ?? ?? ???%模態剛度

omega0=diag(sqrt(Vnorm'*K*Vnorm)) ;? ?? ?? ?? ?? ? % 正則化模態剛度

VnormM=diag(Vnorm'*M*Vnorm);? ?? ?? ?? ?? ?? ?? ???%正則化模態質量

Fnorm=Vnorm'*ft0 ;? ?? ?? ?? ?? ?? ?? ?? ?? ?? ?? ?? ? % 模態力矢量

Modamp=Vnorm'*(a*M+b*K)*Vnorm;? ?? ?? ?? ?? ?? ?? ?% 模態阻尼矩陣

zeta=diag((1/2)*Modamp*inv(diag(omega0)));? ?? ?? ???% 阻尼比

%--------------------------------------------------------------------------

%??%??(4) 模態坐標的響應

%--------------------------------------------------------------------------

q0=zeros(length(K),1);

dq0=zeros(length(K),1);? ?? ?? ?? ?? ?? ?? ?? ? % 初始化位移和速度

eta0=Vnorm'*M*q0; deta0=Vnorm'*M*dq0;? ???%2.2-53% 初始條件模態坐標的位移和速度

eta=zeros(length(t),length(K));? ?? ?? ?? ?? ? %??%初始化位移2.2-52

for i=1:length(K)? ?? ?? ?? ?? ?? ?? ?? ?? ?? ?? ?% t(i)時刻的響應

phase0=omega0(i)*t;? ?? ?? ?? ?? ?? ?? ?? ?? ?? ???%w*t omega0為無阻尼固有頻率

omegad=omega(i)*sqrt(1-zeta(i)^2);

phase=omegad*t;? ?? ?? ?? ?? ?? ?? ?? ?? ? %wd*t其中(omegad為有阻尼固有頻率)

Exx=exp(-zeta(i)*omega(i)*t);? ?? ?? ?? ?? ? %中間變量

C1=eta0(i);? ?? ?? ?? ?? ?? ?? ?? ?? ?? ?? ?? ?? ? %初始位移

C2=(deta0(i)+eta0(i)*zeta(i)*omega0(i))/omegad;? ?? ?? ? %中間變量

D1=zeta(i)*omega(i)/omegad;

II=ones(length(t),1);

XX=Fnorm(i)/(omegad^2+zeta(i)^2*omega(i)^2);

eta(:,i)=C1*Exx.*cos(phase)+C2*Exx.*sin(phase)+ XX*(II-Exx.*cos(phase)-D1*Exx.*sin(phase));%有阻尼系統

%??eta(:,i)=Fnorm(i)/((omega(i))^2)*(1-cos(omega(i)*t)); %無阻尼系統

end

%--------------------------------------------------------------------------

%? ?(6) 將模態坐標轉化到物理坐標系

%--------------------------------------------------------------------------

eta=eta';? ?? ?? ?? ?%模態坐標

y=Vnorm*eta;? ?? ?? ???%物理坐標

plot(t, y(2,:),'-r*')

xlabel('Time(seconds)')

ylabel(' Vertical displ. (m)')

總結

以上是生活随笔為你收集整理的matlab求微分方程的系数,如何利用matlab求解矩阵系数的二阶微分方程的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。