MATLAB用newmarkβ和rungekutta求解duffing振子和非线性加幂律振子响应.zip

%% 生成白噪声样本 clear all; close all; N=10000;%噪声点数 B=100;%带宽100Hz； dt=1/B; t=0:dt:(N-1)*dt; zeta=0.05;%阻尼比为0.05； Sw=2*zeta/pi;%功率谱密度； f=0-B/2:B/N:B-B/N-B/2;%频谱列表 w=wgn(1,N,Sw*B,'linear');%功率P = 功率谱密度*B = k*B； figure(1) plot(t,w);%绘制出一万个样本的图 title('样本函数') %% y1=fft(w,N); p1=(y1.*conj(y1))/N/B; %直接法估计功率谱（对样本进行fft再取模的平方除以点数） figure(2); hist(w,100); title('白噪声的概率分布'); histfit(w); legend('频率分布直方图','拟合概率分布曲线') figure(3) plot(f,p1); title('白噪声功率谱'); xlabel('f/Hz'); ylabel('Sw'); fprintf('噪声功率谱密度均值为：%f\n', mean(p1/N)/B); %% 设定初始值 m=1; wn=1; c=2*zeta*m*wn; k=m*wn^2; m_size=size(m); DOF=m_size(1);%自由度 p=w; Npoint=length(p); DeltT=1/B; epsilon=0.1; [u1,uu1,uuu1] = newmarkduffing(m,c,k,DOF,epsilon,DeltT,Npoint,p); [u2,uu2,uuu2] = newmarkjiamilv(m,c,k,DOF,epsilon,DeltT,Npoint,p); figure(4) plot(t,u1); xlabel('time/s');ylabel('displacement/m'); hold on plot(t,u2); legend('duffing振子','非线性加幂律阻尼振子'); title('Newmark β求解位移') hold off pt=t; [T,Y] = ode45(@(t,y) myode(t,y,pt,p,c,k,epsilon,m),[0 t(end)],[0;0]); figure(5) plot(t,u1,T,Y(:,1),'-') title('Runge-kutta与Newmark β求解duffing振子'); legend('duffing振子','非线性加幂律阻尼振子'); xlabel('time/s');ylabel('displacement/m'); [T1,Y1] = ode45(@(t,y) myode2(t,y,pt,p,c,k,epsilon,m),[0 t(end)],[0;0]); figure(6) plot(t,u2,T1,Y1(:,1),'-'); title('Runge-kutta与Newmark β求解非线性加幂律阻尼振子'); legend('duffing振子','非线性加幂律阻尼振子'); xlabel('time/s');ylabel('displacement/m'); %% 对比epsilon对响应的影响 u1=zeros(6,Npoint);u2=zeros(6,Npoint); for i=1:6 epsilon=0:0.2:1; epsilon=epsilon(i); [u1(i,:),uu1,uuu1] = newmarkduffing(m,c,k,DOF,epsilon,DeltT,Npoint,p); [u2(i,:),uu2,uuu2] = newmarkjiamilv(m,c,k,DOF,epsilon,DeltT,Npoint,p); end figure(7) hold on for i=1:6 plot(t,u1(i,:)) end legend('epsilon=0','epsilon=0.2','epsilon=0.4','epsilon=0.6','epsilon=0.8','epsilon=1'); hold off xlabel('time/s');ylabel('displacement/m'); figure(8) hold on for i=1:6 plot(t,u2(i,:)) end legend('epsilon=0','epsilon=0.2','epsilon=0.4','epsilon=0.6','epsilon=0.8','epsilon=1'); xlabel('time/s');ylabel('displacement/m'); hold off