How to solve multiple DOF mass-spring linear system with attached resonators with ode45?

The system is this:

I have the initial conditions, but would like to know how to solve this system with ode45 or any other solver, because they are coupled equations.

Thanks in advance!

Answer 1

[Alan Stevens]

Something like this perhaps (but use your own data!):

% Initial conditions
u = zeros(length(tspan),length(x));
p = u; p(1,N/2) = 1;
v = u; 
q = u;
ic = [u(1,:), p(1,:), v(1,:), q(1,:)];

% Call ode solver
[t, y] = ode45(@rates, tspan, ic);

% Plot results
u = y(:,1:N);  v = y(:,2*N+1:3*N);
plot(t,u),grid
xlabel('time'), ylabel('u')
figure
plot(t,v),grid
xlabel('time'), ylabel('v')

function dydt = rates(~,y)
        % arbitrary data
        m1 = 5; m2 = 1;
        k1 = 1; k2 = 2;
        
        n = length(y)/4;
        
        u = y(1:n); p = y(n+1:2*n);
        v = y(2*n+1:3*n); q = y(3*n+1:end);
        
        dudt(1) = p(1);
        dpdt(1) = (k1/m1)*(-u(1)+u(2)) + (u(1)-v(1))/m1;
        dvdt(1) = q(1);
        dqdt(1) = (k2/m2)*(u(1)-v(1));
        
        % calculate rates
        for j = 2:n-1

            dudt(j) = p(j);
            dpdt(j) = (k1/m1)*(u(j-1)-2*u(j)+u(j+1)) + (u(j)-v(j))/m1;

            dvdt(j) = q(j);
            dqdt(j) = (k2/m2)*(u(j)-v(j));

        end
         
        dudt(n) = p(n);
        dpdt(n) = (k1/m1)*(-u(n-1)+u(n)) + (u(n)-v(n))/m1;
        dvdt(n) = q(n);
        dqdt(n) = (k2/m2)*(u(n)-v(n));

        dydt = [dudt'; dpdt'; dvdt'; dqdt'];

end

[Dane Stan]

What if I have a prescribed harmonic displacement applied in the middle, i.e. u(n/2)=cos(t)=f(t) (n-odd) where should I write it in the code? I tried

dudt((n+1)/2) = p((n+1)/2);

dpdt((n+1)/2) = (k1/m1)*(u((n+1)/2-1)-2*f(t)+u((n+1)/2+1)) + (f(t)-v((n+1)/2))/m1;

dvdt((n+1)/2) = q((n+1)/2);

dqdt((n+1)/2) = (k2/m2)*(f(t)-v((n+1)/2));

but I think I am not doing it right because I am not getting the desired results.

Thanks again!

[Alan Stevens]

If it's just applied to the u'' equation then perhaps like the following (assuming n is even):

N = 10;             % Number of nodes
x = 0:N-1;          % x-values of nodes
tspan = 0:0.1:5;    % Time values

% Initial conditions
u = zeros(length(tspan),length(x));
p = u; p(1,N/2) = 1;
v = u; 
q = u;
ic = [u(1,:), p(1,:), v(1,:), q(1,:)];

% Call ode solver
[t, y] = ode45(@rates, tspan, ic);

% Plot results
u = y(:,1:N);  v = y(:,2*N+1:3*N);
plot(t,u),grid
xlabel('time'), ylabel('u')
figure
plot(t,v),grid
xlabel('time'), ylabel('v')

function dydt = rates(t,y)  %%%%%%% Put t as an explicit input
        % arbitrary data
        m1 = 5; m2 = 1;
        k1 = 1; k2 = 2;
        
        n = length(y)/4;        
        u = y(1:n); p = y(n+1:2*n);
        v = y(2*n+1:3*n); q = y(3*n+1:end);
        
        dudt(1) = p(1);
        dpdt(1) = (k1/m1)*(-u(1)+u(2)) + (u(1)-v(1))/m1;
        dvdt(1) = q(1);
        dqdt(1) = (k2/m2)*(u(1)-v(1));
        
        % calculate rates
        for j = 2:n-1

            dudt(j) = p(j);
            dpdt(j) = (k1/m1)*(u(j-1)-2*u(j)+u(j+1)) + (u(j)-v(j))/m1;
            dvdt(j) = q(j);
            dqdt(j) = (k2/m2)*(u(j)-v(j));

        end
        dpdt(n/2) = dpdt(n/2) + cos(t)/m1; %%%%% forcing function on midpoint 
         
        dudt(n) = p(n);
        dpdt(n) = (k1/m1)*(-u(n-1)+u(n)) + (u(n)-v(n))/m1;
        dvdt(n) = q(n);
        dqdt(n) = (k2/m2)*(u(n)-v(n));

        dydt = [dudt'; dpdt'; dvdt'; dqdt'];

end

I'd find it easier to decide if you wrote the mathematical equations (rather than the computer ones) including the cos(t) forcing function.