function zdzdxq m = 20; k = 10; L0 = 10; w0 = 0.1; x0 = 1; tmax = 50; h = 0.05; [t, u] = ode45(@F, 0:h:tmax, [x0 0 0 w0]); % 需定义函数F subplot(2,2,1) L = L0 + 2*u(:,1); plot(t, L, 'linewidth', 1.5), title('弹簧长度-时间图像'); xlabel('t'); ylabel('l'); axis([0 tmax 8 12.5]); subplot(2,2,2) theta = u(:,3) - floor(u(:,3)/(2*pi)) * 2*pi; % 修正括号和运算符 plot(t, theta, 'linewidth', 1.5), title('角度-时间图像'); xlabel('t'); ylabel('\theta'); axis([0 tmax 0 2*pi+0.1]); subplot(2,2,3) plot(t, u(:,4), 'linewidth', 1.5), title('角速度-时间图像'); xlabel('t'); ylabel('\omega'); axis([0 tmax 0 0.3]); subplot(2,2,4) x1 = (L0/2 + u(:,1)) .* cos(u(:,3)); y1 = (L0/2 + u(:,1)) .* sin(u(:,3)); x2 = (L0/2 - u(:,1)) .* cos(u(:,3)); y2 = (L0/2 - u(:,1)) .* sin(u(:,3)); xx = linspace(x1(1), x2(1), 50); % 弹簧两端应为x1和x2 yy = 0.35*sin(20*pi*(xx - x2(1))/L0) + y1(1); xx1 = (xx - x2(1)) .* cos(u(1,3)) - (yy - y1(1)) .* sin(u(1,3)) + x2(1); yy1 = (xx - x2(1)) .* sin(u(1,3)) + (yy - y1(1)) .* cos(u(1,3)) + y1(1); hold on h1 = plot(x1(1), y1(1), 'r.', 'MarkerSize', 45); h2 = plot(x2(1), y2(1), 'g.', 'MarkerSize', 45); % 修正y2(1) h3 = plot(xx1, yy1, 'b', 'LineWidth', 2); axis equal; axis([-6 6 -6 6]); axis off; title('双振子随时间变化'); for i = 1:length(x1) % 修正循环变量和长度 xx = linspace(x1(i), x2(i), 50); yy = 0.35*sin(20*pi*(xx - xx(i))/L(1)) + y1(i); xx_rot = (xx - x2(i)) .* cos(u(i,3)) - (yy - y1(i)) .* sin(u(i,3)) + x2(i); % 使用u(i,3) yy_rot = (xx - x2(i)) .* sin(u(i,3)) + (yy - y1(i)) .* cos(u(i,3)) + y1(i); set(h1, 'XData', x1(i), 'YData', y1(i)); set(h2, 'XData', x2(i), 'YData', y2(i)); set(h3, 'XData', xx_rot, 'YData', yy_rot); drawnow; end end % 必须定义微分方程函数F function dudt = F(t, u) m = 20; k = 10; L0 = 10; x = u(1); dxdt = u(2); theta = u(3); dthetadt = u(4); % 根据系统动力学编写方程 d2xdt2 = (k/m)*(L0/2 - x) - (k/m)*x - (dxdt/m); % 示例，需根据实际模型修正 d2thetadt2 = - (2*x*dthetadt + k*theta)/m; % 示例，需根据实际模型修正 dudt = [dxdt; d2xdt2; dthetadt; d2thetadt2]; end