拉普拉斯变换解mx''+cx'+kx=|sin(t)|遇MATLAB求解问题咨询

Hey there! Let's tackle your problem with the differential equation and MATLAB solving hurdles you're facing.

Is dsolve suitable for this equation? #

Short answer: Yes, it is, but you need to account for the piecewise, periodic nature of |sin(t)| when setting up the problem. The default dsolve call might struggle because it doesn't automatically parse the absolute value of a trigonometric function into a form it can handle cleanly. You'll need to explicitly define the piecewise behavior or leverage the periodicity of |sin(t)| (it has a period of π) to help MATLAB's symbolic engine process the equation correctly.

Feasible solutions to solve and plot the response #

Here are three reliable approaches to get a plotable solution:

1. Fix the symbolic solve with dsolve #

Define the piecewise form of |sin(t)| explicitly for MATLAB's symbolic engine. Since |sin(t)| = sin(t) for 2nπ ≤ t < (2n+1)π and -sin(t) for (2n+1)π ≤ t < (2n+2)π (where n is a non-negative integer), you can use piecewise to formalize this:

syms x(t) m c k
% Define piecewise |sin(t)| with periodic extension
f = piecewise( ...
    t >= 0 & t < pi, sin(t), ...
    t >= pi & t < 2*pi, -sin(t), ...
    t >= 2*pi, piecewise( ...
        t - 2*pi >= 0 & t - 2*pi < pi, sin(t - 2*pi), ...
        t - 2*pi >= pi & t - 2*pi < 2*pi, -sin(t - 2*pi) ...
    ) ...
);

% Define the differential equation
eq = m*diff(x,t,2) + c*diff(x,t) + k*x == f;

% Set initial conditions (adjust these if you have non-zero initial state)
conds = [x(0) == 0, diff(x,t,0) == 0];

% Solve symbolically
sol = dsolve(eq, conds);

% Substitute your actual m, c, k values (example values here)
sol_subs = subs(sol, {m, c, k}, {1, 0.5, 2});

% Plot the solution
fplot(sol_subs, [0, 4*pi]);
xlabel('t'); ylabel('x(t)'); title('Symbolic Solution to mẍ + cẋ + kx = |sin(t)|');
grid on;

Note: For longer time spans, the symbolic expression might get complex, but this works well for verifying short-term behavior.

2. Numerical solving with ode45 (most reliable for plotting) #

Symbolic methods can get bogged down with piecewise periodic functions, so numerical ODE solving is often the easiest path to a plotable result. Convert the second-order ODE into a first-order system, then use ode45:

% Define your constants
m = 1; c = 0.5; k = 2;

% Convert second-order ODE to first-order system:
% Let y1 = x, y2 = dx/dt
% Then dy1/dt = y2, dy2/dt = (|sin(t)| - c*y2 - k*y1)/m
dfun = @(t, y) [y(2); (abs(sin(t)) - c*y(2) - k*y(1))/m];

% Initial conditions (x(0)=0, dx/dt(0)=0; adjust as needed)
y0 = [0; 0];

% Time span to solve over
tspan = [0, 4*pi];

% Run numerical solver
[t, y] = ode45(dfun, tspan, y0);

% Plot the displacement x(t)
plot(t, y(:,1));
xlabel('t'); ylabel('x(t)'); title('Numerical Solution to mẍ + cẋ + kx = |sin(t)|');
grid on;

This method will always produce a smooth, plotable curve and is far faster for large time spans.

3. Use Fourier series expansion of |sin(t)| #

Since |sin(t)| is a periodic function, you can expand it into a Fourier series, compute the steady-state response for each harmonic, add the transient response of the homogeneous equation, and sum them up. This gives you a clear breakdown of the solution's components:

m = 1; c = 0.5; k = 2;
t = linspace(0, 4*pi, 1000);

% Transient response (zero initial conditions here)
lambda = roots([m c k]);
transient = exp(lambda(1)*t)*( (lambda(2)*0 - 0)/(lambda(2)-lambda(1)) ) + ...
            exp(lambda(2)*t)*( (lambda(1)*0 - 0)/(lambda(1)-lambda(2)) );

% Steady-state response via Fourier series (sum first N harmonics)
N = 20;  % Adjust for accuracy vs. speed
steady = 2/pi * ones(size(t));
for n = 1:N
    omega = 2*n;
    % Magnitude and phase of steady-state response for each harmonic
    mag = 1/sqrt( (k - m*omega^2)^2 + (c*omega)^2 );
    phase = -atan2(c*omega, k - m*omega^2);
    steady = steady - (4/(pi*(4*n^2 - 1))) * mag * cos(omega*t + phase);
end

% Total solution
x_total = transient + steady;

% Plot
plot(t, x_total);
xlabel('t'); ylabel('x(t)'); title('Solution via Fourier Series Expansion');
grid on;

This approach is great for understanding how each frequency component contributes to the overall response.

内容的提问来源于stack exchange，提问作者W. Ljungbergh