在matlab中用龙格库塔法求解一个强制质量-弹簧-阻尼系统答案

【问题标题】：Solving a forced mass-spring-damper system with Runge Kutta method in matlab在matlab中用龙格库塔法求解一个强制质量-弹簧-阻尼系统
【发布时间】：2014-02-05 16:06:16
【问题描述】：

我正在尝试使用 Runge-Kutta 方法在 matlab 中解决强制质量弹簧阻尼系统。目前，代码使用常量值作为系统输入，但我想将向量作为输入。例如，我想用矢量值替换我的力幅值 F0。

我应该使用 for 循环还是最简单的方法？

function O = MSDSRK(m,b,k,F0,w,x0,v0) 

% ----- Input argument -----
% m: mass for particle 
% b: damping coefficient 
% k: spring constant 
% F0: amplitude of external force 
% w: angular freuency of external force 
% x0: initial condition for the position x(0)
% v0: initial condition for the velocity v(0)

dt=0.1;

options=odeset('InitialStep',dt,'MaxStep',dt);

td=[0:dt:50];

% Solve differential equation with Runge-Kutta solver 

[t,x]=ode45(@(t,X)MSD(t,X,m,b,k,F0,w),td,[x0;v0],options);

% Extract only particle position trajectory

O=[t x(:,1)];

end


function dX=MSD(t,X,m,b,k,F0,w)

% With two 1st order diffeential equations,
% obtain the derivative of each variables
% (position and velocity of the particle)

dX(1,1)=X(2,1); 

dX(2,1)=(1/m)*(F0*sin(w*t)-b*X(2,1)-k*X(1,1)); 
end

【问题讨论】：

在这种情况下我会使用循环，因为你的函数的输出已经是一个向量。如果你真的想向量化它，让你函数以矩阵的形式产生输出，其中一个维度对应于不同的 F0 值。

标签： matlab

【解决方案1】：

for-loop 肯定会起作用，但在这种情况下我不建议这样做。这种方法的缺点是让你以某种方式思考，也就是说，你正在解决一个简单的一维系统，只是多次。

for-loop 方法所教给您的不仅仅是您已经知道如何做的事情。你会得到一个及格的分数，我敢肯定。但如果你想出类拔萃，不仅在这里，而且在以后的课程中（更重要的是，在以后的工作中），你必须做得更多。您只能通过以不同的方式思考这个简单的问题来学习如何解决更复杂的问题，即作为一个多维质量/弹簧/阻尼器系统，其中所有组件碰巧解耦并且所有发生具有相同的初始值。

对于这样的系统，向量化和矩阵操作确实是要走的路。这是 1D/2D/3D 系统到任意 D 系统的推广。而矩阵操作和向量化，正是 MATLAB 最擅长的。这确实是你可以在这里学到的。

以下是针对您的案例的操作方法。我把它的范围缩小了一点，但原则保持不变：

function O = MSDSRK

    % Assign some random scalar values to all parameters, but 
    % a *vector* to the force amplitudes 
    [m,b,k,F0,w,x0,v0] = deal(...
        1,0.2,3,  rand(13,1),  2,0,0);

    % We need to simulate as many mass-spring-damper systems as there are
    % force amplites, so replicate the initial values
    x0 = repmat(x0, numel(F0),1);
    v0 = repmat(v0, numel(F0),1);

    % The rest is the same as before
    dt = 0.1;
    options = odeset('InitialStep', dt,'MaxStep', dt);

    td = 0:dt:50;

    [t,x] = ode45(@(t,X) MSD(t,X,m,b,k,F0,w), td, [x0;v0], options);

    % The output is now: 
    %
    % x(:,1)  position of mass in system with forcing term F0(1)
    % x(:,2)  position of mass in system with forcing term F0(2)
    % ...
    % x(:,14) speed of mass in system  with forcing term F0(1)
    % x(:,15) speed of mass in system  with forcing term F0(2)
    % ...    
    O = [t x(:,1:numel(F0))];

    % And, to make the differences more clear:
    plot(t, x(:,1:numel(F0)))

end


function dX = MSD(t,X,m,b,k,F0,w)

    % Split input up in position and velocity components
    x = X(1:numel(X)/2);
    v = X(numel(X)/2+1:end);

    % Compute derivative
    dX = [
        v
        (F0*sin(w*t) - b*v - k*x)/m
    ];
end

【讨论】：