求助：用dsolve与ode45复现弹簧阻尼系统Simulink仿真曲线

Let’s get your simulations sorted out—first fixing the dsolve failure you ran into, then breaking down ode45 step by step so it doesn’t feel confusing anymore. We’ll use your provided parameters (m1=500, c1=1200, k1=25000, k2=15000, m2=50) and the step input (step time=2, final value=0.5) to match your Simulink results.

1. Fixing dsolve for Your Two-DOF System #

The first issue in your code is a small syntax typo (eqqn33 vs eqn33) that’s causing dsolve to throw an error. We also need to properly define the step input y(t) instead of leaving it as an undefined symbolic variable. Let’s rebuild this correctly:

Step-by-Step dsolve Code #

syms x(t)
% Define the step input: 0 for t<2, 0.5 for t>=2
y(t) = 0.5 * heaviside(t - 2);
Dy = diff(y, t); % Derivative of step is impulse (heaviside handles this)
Dx = diff(x, t);
D2x = diff(x, 2, t);

% Initial conditions (your original x(0)=0, x'(0)=5—y's conditions are fixed by the input)
cond = [x(0) == 0, Dx(0) == 5];

% Derive the differential equation from your parameters (instead of hardcoded constants)
% From free-body analysis for m1: m1*x'' + c1*x' + (k1+k2)*x = k2*y + c1*y'
% Normalize by m1 to match your original equation structure:
c1_over_m1 = c1/m1;
k_sum_over_m1 = (k1 + k2)/m1;
k2_over_m1 = k2/m1;

eqn = D2x + c1_over_m1*Dx + k_sum_over_m1*x == k2_over_m1*y + c1_over_m1*Dy;

% Solve the equation
sol = dsolve(eqn, cond);

% Simplify and plot the result
x_sol = simplify(sol);
fplot(x_sol, [0, 10]); % Plot from t=0 to 10 seconds
xlabel('Time (s)');
ylabel('Displacement x(t)');
title('dsolve Solution for Spring-Damper System');
grid on;

This should resolve your dsolve failure—we fixed the typo, defined the input properly, and used your actual parameters to avoid hardcoded constants that might not align with your system.

2. Understanding ode45 for Your System #

ode45 solves systems of first-order ordinary differential equations, so we need to convert your second-order equation into a set of first-order equations. Here’s how to do it clearly:

Step 1: Convert to First-Order State Variables #

Let’s define two state variables to represent our system:

• z1 = x (displacement of mass m1)
• z2 = x' (velocity of mass m1)

Now, we can write the derivatives of these states:

• z1' = z2 (derivative of displacement is velocity)
• z2' = (-c1*z2 - (k1+k2)*z1 + k2*y(t) + c1*y'(t))/m1 (from rearranging the second-order equation)

Step 2: Write the ODE Function #

Create a separate function file (name it spring_damper_ode.m) that computes the state derivatives:

function dzdt = spring_damper_ode(t, z, m1, c1, k1, k2, y_final)
    z1 = z(1); % Displacement x
    z2 = z(2); % Velocity x'
    
    % Handle the step input
    if t < 2
        y = 0;
        Dy = 0;
    else
        y = y_final;
        Dy = 0; % After t=2, the step input has zero derivative
    end
    
    % Calculate the derivatives
    dzdt = [z2;
            (-c1*z2 - (k1 + k2)*z1 + k2*y + c1*Dy)/m1];
end

Step 3: Run the ode45 Simulation #

In your main script, we’ll split the simulation into two parts (before and after the step at t=2) to avoid dealing with the impulse from the step’s derivative:

% System parameters
m1 = 500;
c1 = 1200;
k1 = 25000;
k2 = 15000;
y_final = 0.5;

% Initial conditions: x(0)=0, x'(0)=5
z0 = [0; 5];

% First interval: t=0 to 2 (input is 0)
t_interval1 = [0, 2];
[ts1, zs1] = ode45(@(t,z) spring_damper_ode(t,z,m1,c1,k1,k2,y_final), t_interval1, z0);

% Second interval: t=2 to 10 (input is 0.5), use final state from first interval as initial
z0_interval2 = zs1(end,:)';
t_interval2 = [2, 10];
[ts2, zs2] = ode45(@(t,z) spring_damper_ode(t,z,m1,c1,k1,k2,y_final), t_interval2, z0_interval2);

% Combine results (avoid duplicate t=2)
t_total = [ts1; ts2(2:end)];
x_total = [zs1(:,1); zs2(2:end,1)];

% Plot the result
plot(t_total, x_total);
xlabel('Time (s)');
ylabel('Displacement x(t)');
title('ode45 Solution for Spring-Damper System');
grid on;

Why This Works #

By splitting the simulation at t=2, we avoid the messy impulse from the step’s derivative, which ode45 handles cleanly. The first part uses the initial input (0), and the second part uses the final input value (0.5) with no derivative change.

3. Matching Simulink Results #

To ensure your Simulink model aligns with these results:

• Double-check that your Simulink blocks use exactly the same parameters (m1=500, c1=1200, etc.)
• Confirm your step input block is set to Step time=2, Final value=0.5, Initial value=0
• Verify initial conditions in Simulink: set the integrator for x to 0, and the integrator for x’ to 5

If you’re working with the second system (I assume it’s a different configuration, maybe single-DOF?), follow the same process: derive its differential equation, convert to first-order states for ode45, or use dsolve with the correct input definition.

内容的提问来源于stack exchange，提问作者S. Seo