Решение системы принудительной массы-пружины-демпфера методом Рунге-Кутты в Matlab

Я пытаюсь решить систему принудительного масс-пружинного демпфера в 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

matlab
person user3200392    schedule 04.02.2014    source источник
comment
В этом случае я бы использовал цикл, поскольку вывод вашей функции уже является вектором. Если вы действительно хотите его векторизовать, заставьте свою функцию выводить результат в виде матрицы, одно из измерений которой соответствует разным значениям F0.   -  person freude    schedule 04.02.2014

Ответы (1)


arrow_upward
0
arrow_downward

Цикл for точно сработает, но в данной ситуации я бы не советовал его использовать. Недостаток этого подхода состоит в том, что он заставляет вас думать определенным образом, а именно, что вы решаете простую одномерную систему всего несколько раз.

Подход с циклом for не научит вас большему, чем вы уже знаете, как это делать. Ты получишь проходной балл, я уверен. Но если вы хотите преуспеть не только здесь, но и на будущих занятиях (и, что более важно, на работе позже), вы должны делать больше. Вы можете научиться решать более сложные проблемы, только рассматривая эту простую проблему по-другому, а именно как многомерную систему масса/пружина/демпфер, все компоненты которой случайно являются несвязаны, и все они случайно имеют одинаковые начальные значения.

Для таких систем векторизация и манипулирование матрицами — правильный путь. Это обобщение систем 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
person Rody Oldenhuis    schedule 05.02.2014