matlab解弹簧谐振子二阶齐次微分方程

matlab解弹簧谐振子二阶齐次微分方程spring-oscillator-ode-solverMATLAB解弹簧谐振子二阶齐次微分方程V1matlab%弹簧谐振子二阶齐次微分方程求解clearall;closeall;clc;%定义系统参数m=1;%质量(kg)k=10;%弹簧刚度(N/m)c=0.5;%阻尼系数(N·s/m)x0=1;%初始位置(m)v0=0;%初始速度(m/s)tSpan=[010];%时间范围(s)%定义微分方程dydt=@(t,y)[y(2);%y(1)的导数是y(2)（速度）-(k/m)*y(1)-(c/m)*y(2)];%y(2)的导数是加速度%设置初始条件y0=[x0;v0];%使用ode45求解[t,y]=ode45(dydt,tSpan,y0);%计算理论解（无阻尼情况对比）wn=sqrt(k/m);%自然频率ifc==0%无阻尼情况x_theory=x0*cos(wn*t)+(v0/wn)*sin(wn*t);else%有阻尼情况，根据阻尼比计算zeta=c/(2*sqrt(m*k));%阻尼比ifzeta<1%欠阻尼wd=wn*sqrt(1-zeta^2);%阻尼自然频率x_theory=exp(-zeta*wn*t).*(x0*cos(wd*t)+((v0+zeta*wn*x0)/wd)*sin(wd*t));elseifzeta==1%临界阻尼x_theory=(x0+(v0+wn*x0)*t).*exp(-wn*t);else%过阻尼s1=-wn*(zeta+sqrt(zeta^2-1));s2=-wn*(zeta-sqrt(zeta^2-1));A=(v0-s2*x0)/(s1-s2);B=x0-A;x_theory=A*exp(s1*t)+B*exp(s2*t);endend%绘制结果figure('Position',[100,100,1000,400])%位置-时间图subplot(1,2,1)plot(t,y(:,1),'b-','LineWidth',1.5)holdonplot(t,x_theory,'r--','LineWidth',1.5)xlabel('时间(s)')ylabel('位置(m)')title('弹簧谐振子位移响应')legend('数值解','理论解')gridon%相图subplot(1,2,2)plot(y(:,1),y(:,2),'g-','LineWidth',1.5)xlabel('位置(m)')ylabel('速度(m/s)')title('弹簧谐振子相图')gridon%打印系统参数和结果fprintf('系统参数:\n')fprintf('质量:%.2fkg\n',m)fprintf('弹簧刚度:%.2fN/m\n',k)fprintf('阻尼系数:%.2fN·s/m\n',c)fprintf('阻尼比:%.4f\n',c/(2*sqrt(m*k)))fprintf('自然频率:%.4frad/s\n',sqrt(k/m))这个程序解决了一般形式的弹簧-质量-阻尼系统的运动方程：mx¨+cx˙+kx=0其中：m

是质量c

是阻尼系数k

是弹簧刚度程序将二阶微分方程转换为一阶方程组，然后使用MATLAB的ode45求解