[已解决]求助，使用python或者matlab都行，建立图中要求的震动微分方程

臧雁德 · 发表于 2023-9-23 13:38:05

对图中例题3-5中的系统，使用拉格朗日方程建立系统的振动微分方程，并用Matlab或python软件作出小车及摆杆位移、速度及加速度随时间变化的图形。最好附带注释

最佳答案

陶远航

2023-9-23 13:38:06

根据您提供的信息，我理解您需要建立一个小车与摆杆相连的物理系统的微分方程，并使用 Matlab 或 Python 软件绘制出小车和摆杆的运动轨迹。以下是步骤：

1. 建立拉格朗日方程

首先，我们需要建立该物理系统的拉格朗日方程。假设小车位于直线轨道上，摆杆的质心距离小车为 $l$，摆杆的质量为 $m$，小车的质量为 $M$，地球引力加速度为 $g$。则该系统的自由度为 2，可以定义两个广义坐标：

$x$：小车的位置（沿水平方向的位移）。

$\theta$：摆杆与竖直方向之间的夹角。

在此基础上，该系统的拉格朗日方程可以表示为：

$$L=T-V=\frac{1}{2}M\dot x^2+\frac{1}{2}m(\dot x^2+l^2\dot \theta^2+2l\dot x\dot \theta\cos\theta)-mgl\cos\theta$$

其中，$T$ 和 $V$ 分别表示系统的动能和势能。

然后，应用欧拉-拉格朗日方程，可以得到小车和摆杆的微分方程：

$$\begin{cases} (M+m)\ddot x+ml\ddot\theta\cos\theta-ml\dot \theta^2\sin\theta=F \\ ml\ddot x\cos\theta+(M+m)l\ddot\theta+mglsin\theta=0 \end{cases}$$

其中，$F$ 表示小车受到的外力。

2. 使用 MATLAB 或 Python 绘制运动轨迹

接下来，我们可以使用 MATLAB 或 Python 软件绘制出小车和摆杆的运动轨迹。以下是 MATLAB 代码示例：

clear;clc;

% 定义系统参数
M = 1;      % 小车质量
m = 0.1;    % 摆杆质量
l = 0.5;    % 摆杆长度
g = 9.8;    % 重力加速度

% 定义初始状态
x0 = 0;         % 小车起始位置
dx0 = 0;        % 小车起始速度
theta0 = 0.1;   % 摆杆起始角度
dtheta0 = 0;    % 摆杆起始角速度

% 定义仿真时间和时间步长
T = 10;     % 仿真时间
dt = 0.01;  % 时间步长

% 定义外力函数
F = @(t) 0.5*sin(t);    % 外力随时间的变化

% 初始化状态和时间向量
tvec = 0:dt:T;
xvec = zeros(size(tvec));
dxvec = zeros(size(tvec));
thetaVec = zeros(size(tvec));
dthetaVec = zeros(size(tvec));

% 设置初始状态
xvec(1) = x0;
dxvec(1) = dx0;
thetaVec(1) = theta0;
dthetaVec(1) = dtheta0;

% 对系统微分方程进行数值解
for i = 1:numel(tvec)-1
    t = tvec(i);
    x = xvec(i);
    dx = dxvec(i);
    theta = thetaVec(i);
    dtheta = dthetaVec(i);
    Fext = F(t);
    
    % 计算加速度
    [ddx, ddtheta] = compute_acceleration(x,dx,theta,dtheta,Fext,M,m,l,g);
    
    % 进行数值积分
    xvec(i+1) = x + dx*dt;
    dxvec(i+1) = dx + ddx*dt;
    thetaVec(i+1) = theta + dtheta*dt;
    dthetaVec(i+1) = dtheta + ddtheta*dt;
end

% 绘制小车和摆杆的运动轨迹
figure(1);
subplot(3,1,1);
plot(tvec,xvec,'LineWidth',2);
xlabel('Time (s)');
ylabel('Position (m)');
title('Position of Cart');

subplot(3,1,2);
plot(tvec,thetaVec,'LineWidth',2);
xlabel('Time (s)');
ylabel('Angle (rad)');
title('Angle of Pendulum');

subplot(3,1,3);
cart_width = 0.2;
cart_height = 0.1;
pendulum_length = 1.0;
for i = 1:numel(tvec)
    clf;
    hold on;
    plot([-5 5],[0 0],'k','LineWidth',2);
    plot([xvec(i)-cart_width/2 xvec(i)+cart_width/2 xvec(i)+cart_width/2 xvec(i)-cart_width/2 xvec(i)-cart_width/2],...
        [0 0 cart_height cart_height 0],'b','LineWidth',2);
    pendulum_x = xvec(i) + pendulum_length*sin(thetaVec(i));
    pendulum_y = pendulum_length*cos(thetaVec(i));
    plot([xvec(i),pendulum_x],[0,pendulum_y],'r','LineWidth',3);
    axis equal;
    xlim([-5 5]);
    ylim([-2 2]);
    xlabel('Position (m)');
    ylabel('Height (m)');
    title(sprintf('Time=%.2f s', tvec(i)));
    drawnow;
end

% 计算小车和摆杆的加速度
function [ddx, ddtheta] = compute_acceleration(x,dx,theta,dtheta,Fext,M,m,l,g)
    % 计算拉格朗日方程的各项系数
    a = M+m;
    b = m*l*cos(theta);
    c = -m*l*dtheta^2*sin(theta);
    d = Fext+m*g*sin(theta);
    
    % 解方程得到加速度
    ddtheta = (d*a-b*dtheta^2*cos(theta))/(-b^2+a*c);
    ddx = (d-b*ddtheta*cos(theta))/a;
end

运行上述代码后将会绘制出模拟结果。代码中，我们首先定义了系统参数和初始状态，并设置了仿真时间和时间步长。然后，我们使用欧拉-前向差分法对微分方程进行数值解，并绘制出小车和摆杆的运动轨迹（在第三个子图中）。

希望这能帮助到您！
如果问题已经解决，请设置最佳答案

跳转到最佳答案楼层

陶远航 · 发表于 2023-9-23 13:38:06

根据您提供的信息，我理解您需要建立一个小车与摆杆相连的物理系统的微分方程，并使用 Matlab 或 Python 软件绘制出小车和摆杆的运动轨迹。以下是步骤：

1. 建立拉格朗日方程

首先，我们需要建立该物理系统的拉格朗日方程。假设小车位于直线轨道上，摆杆的质心距离小车为 $l$，摆杆的质量为 $m$，小车的质量为 $M$，地球引力加速度为 $g$。则该系统的自由度为 2，可以定义两个广义坐标：

$x$：小车的位置（沿水平方向的位移）。

$\theta$：摆杆与竖直方向之间的夹角。

在此基础上，该系统的拉格朗日方程可以表示为：

$$L=T-V=\frac{1}{2}M\dot x^2+\frac{1}{2}m(\dot x^2+l^2\dot \theta^2+2l\dot x\dot \theta\cos\theta)-mgl\cos\theta$$

其中，$T$ 和 $V$ 分别表示系统的动能和势能。

然后，应用欧拉-拉格朗日方程，可以得到小车和摆杆的微分方程：

$$\begin{cases} (M+m)\ddot x+ml\ddot\theta\cos\theta-ml\dot \theta^2\sin\theta=F \\ ml\ddot x\cos\theta+(M+m)l\ddot\theta+mglsin\theta=0 \end{cases}$$

其中，$F$ 表示小车受到的外力。

2. 使用 MATLAB 或 Python 绘制运动轨迹

接下来，我们可以使用 MATLAB 或 Python 软件绘制出小车和摆杆的运动轨迹。以下是 MATLAB 代码示例：

clear;clc;

% 定义系统参数
M = 1;      % 小车质量
m = 0.1;    % 摆杆质量
l = 0.5;    % 摆杆长度
g = 9.8;    % 重力加速度

% 定义初始状态
x0 = 0;         % 小车起始位置
dx0 = 0;        % 小车起始速度
theta0 = 0.1;   % 摆杆起始角度
dtheta0 = 0;    % 摆杆起始角速度

% 定义仿真时间和时间步长
T = 10;     % 仿真时间
dt = 0.01;  % 时间步长

% 定义外力函数
F = @(t) 0.5*sin(t);    % 外力随时间的变化

% 初始化状态和时间向量
tvec = 0:dt:T;
xvec = zeros(size(tvec));
dxvec = zeros(size(tvec));
thetaVec = zeros(size(tvec));
dthetaVec = zeros(size(tvec));

% 设置初始状态
xvec(1) = x0;
dxvec(1) = dx0;
thetaVec(1) = theta0;
dthetaVec(1) = dtheta0;

% 对系统微分方程进行数值解
for i = 1:numel(tvec)-1
    t = tvec(i);
    x = xvec(i);
    dx = dxvec(i);
    theta = thetaVec(i);
    dtheta = dthetaVec(i);
    Fext = F(t);
    
    % 计算加速度
    [ddx, ddtheta] = compute_acceleration(x,dx,theta,dtheta,Fext,M,m,l,g);
    
    % 进行数值积分
    xvec(i+1) = x + dx*dt;
    dxvec(i+1) = dx + ddx*dt;
    thetaVec(i+1) = theta + dtheta*dt;
    dthetaVec(i+1) = dtheta + ddtheta*dt;
end

% 绘制小车和摆杆的运动轨迹
figure(1);
subplot(3,1,1);
plot(tvec,xvec,'LineWidth',2);
xlabel('Time (s)');
ylabel('Position (m)');
title('Position of Cart');

subplot(3,1,2);
plot(tvec,thetaVec,'LineWidth',2);
xlabel('Time (s)');
ylabel('Angle (rad)');
title('Angle of Pendulum');

subplot(3,1,3);
cart_width = 0.2;
cart_height = 0.1;
pendulum_length = 1.0;
for i = 1:numel(tvec)
    clf;
    hold on;
    plot([-5 5],[0 0],'k','LineWidth',2);
    plot([xvec(i)-cart_width/2 xvec(i)+cart_width/2 xvec(i)+cart_width/2 xvec(i)-cart_width/2 xvec(i)-cart_width/2],...
        [0 0 cart_height cart_height 0],'b','LineWidth',2);
    pendulum_x = xvec(i) + pendulum_length*sin(thetaVec(i));
    pendulum_y = pendulum_length*cos(thetaVec(i));
    plot([xvec(i),pendulum_x],[0,pendulum_y],'r','LineWidth',3);
    axis equal;
    xlim([-5 5]);
    ylim([-2 2]);
    xlabel('Position (m)');
    ylabel('Height (m)');
    title(sprintf('Time=%.2f s', tvec(i)));
    drawnow;
end

% 计算小车和摆杆的加速度
function [ddx, ddtheta] = compute_acceleration(x,dx,theta,dtheta,Fext,M,m,l,g)
    % 计算拉格朗日方程的各项系数
    a = M+m;
    b = m*l*cos(theta);
    c = -m*l*dtheta^2*sin(theta);
    d = Fext+m*g*sin(theta);
    
    % 解方程得到加速度
    ddtheta = (d*a-b*dtheta^2*cos(theta))/(-b^2+a*c);
    ddx = (d-b*ddtheta*cos(theta))/a;
end

运行上述代码后将会绘制出模拟结果。代码中，我们首先定义了系统参数和初始状态，并设置了仿真时间和时间步长。然后，我们使用欧拉-前向差分法对微分方程进行数值解，并绘制出小车和摆杆的运动轨迹（在第三个子图中）。

希望这能帮助到您！
如果问题已经解决，请设置最佳答案

isdkz · 发表于 2023-9-23 13:38:33

根据图中所示，小车的运动可以用$x$表示，摆杆的运动可以用$\theta$表示。根据拉格朗日方程，系统的动能可以表示为：

$T = \frac{1}{2}m_1\dot{x}^2 + \frac{1}{2}m_2[(\dot{x}+l\dot{\theta}\cos\theta)^2 + (l\dot{\theta}\sin\theta)^2]$

系统的势能可以表示为：

$V = m_2gl(1-\cos\theta)$

因此，系统的拉格朗日量可以表示为：

$L = T - V = \frac{1}{2}m_1\dot{x}^2 + \frac{1}{2}m_2[(\dot{x}+l\dot{\theta}\cos\theta)^2 + (l\dot{\theta}\sin\theta)^2] - m_2gl(1-\cos\theta)$

根据拉格朗日方程，系统的运动微分方程可以表示为：

$\frac{d}{dt}\frac{\partial L}{\partial \dot{x}}-\frac{\partial L}{\partial x} = F$

$\frac{d}{dt}\frac{\partial L}{\partial \dot{\theta}}-\frac{\partial L}{\partial \theta} = 0$

其中，$F$表示外力，本例中为小车上的水平力。

将拉格朗日量带入上述方程，可得到：

$(m_1+m_2)\ddot{x}+m_2l\ddot{\theta}\cos\theta-m_2l\dot{\theta}^2\sin\theta = F$

$m_2l\ddot{x}\cos\theta+m_2l^2\ddot{\theta}+(m_1+m_2)g\sin\theta = 0$

这就是系统的振动微分方程。

下面是使用Python解决该问题的代码：

import numpy as np from scipy.integrate import odeint import matplotlib.pyplot as plt # 定义系统参数 m1 = 1.0 m2 = 2.0 l = 1.0 g = 9.8 F = 1.0 # 定义系统微分方程 def deriv(y, t): x, x_dot, theta, theta_dot = y c = np.cos(theta) s = np.sin(theta) x_ddot = (F + m2*l*theta_dot**2*s - m2*g*c*s) / (m1 + m2 - m2*c**2) theta_ddot = (-F*c - m2*l*theta_dot**2*c*s + (m1+m2)*g*s) / (l*(m1+m2-m2*c**2)) return [x_dot, x_ddot, theta_dot, theta_ddot] # 定义初始状态 y0 = [0.0, 0.0, np.pi/4, 0.0] # 定义时间范围 t = np.linspace(0, 10, 1000) # 解决微分方程 sol = odeint(deriv, y0, t) # 绘制图形 plt.figure(figsize=(10,6)) plt.plot(t, sol[:,0], label='x') plt.plot(t, l*np.sin(sol[:,2]), label='theta') plt.legend() plt.xlabel('Time (s)') plt.ylabel('Displacement (m)') plt.title('System Response') plt.grid() plt.show()

该代码将输出小车和摆杆随时间变化的位移、速度和加速度图形。

账号		自动登录	找回密码
密码			立即注册