Change Input Force to Input Displacement ODE Solver

I have the following ODE solution:

N       = 100;
m       = 0.1*ones(N,1); 
c       = 0.1; 
b       = 0.1; 
k       = 4; 
gamma   = 0.1; 
X0   = zeros(2*N, 1);
dt = 0.91; % [s] 
scale = 0.0049/2; 
epsilon = 0.5; % [m] 
escale = 10^-2;
rd = 1;
f =  @(t,rd) rd.*scale.*square(t) + epsilon.*escale; % [N]
fun = @(t, X) odefun(t, rd, X, N, marray(m), make_diagonal(X,c,b,N),...
make_diagonal(X, k, gamma, N),f);
tspan_train   = [0:dt:100];
[t, X]  = ode45(fun, tspan_train, X0);
function dX = odefun(t, rd, X, N, M, C, K, f)
    %% Definitions
    x       = X(1:N);               % position state vector
    dx      = X(N+1:2*N);           % velocity state vector
    %% Force vector
    f_reservoir       = f(t,rd);
    F       = [f_reservoir; zeros(98,1); f_reservoir];
    %% Equations of Motion
    ddx     = M\(F - K*x - C*dx);
    
    %% State-space model
    dX      = [dx; ...
               ddx];
end
function out = make_diagonal(x, k, gamma, N)
    x = x(:);
    x = [x(1:N); 0];
    ck = circshift(k, -1);
    cg = circshift(gamma, -1);
    cx = circshift(x, -1);
    ccx = circshift(x, 1);
    d1 = -3 .* ck .* cg .* cx .^ 2 - ck;
    d2 = (k .* gamma + ck .* cg) .* x .^ 2 + k + ck;
    d3 = -3 .* k .* ccx .^ 2 - k;
    out = full(spdiags([d1 d2 d3], -1:1, N, N));
end
function M = marray(m)
    M = diag(m);
end

I would like to change my input force f in Newtons to an input displacement in meters. How do I do this?

[Aquatris]

From what I understand you have a system of 100 mass-spring-dampers connected in series. So when you say you want to change your input from force to displacement, where are you planning to apply this displacement input, which mass?

[Jonathan Frutschy]

@Aquatris Yes. These are 100 sprig mass dampers in series. I would like to apply the input displacement to masses 1&100 (end masses). Currently, I apply f or f_reservoir after substitution to masses 1&100.

[Aquatris]

I think you can simply do

function dX = odefun(t, rd, X, N, M, C, K, f,displacementInput)
%% Definitions
x       = X(1:N);               % position state vector
x(1) = x(1)+displacementInput(1); 
x(end) = x(end)+displacementInput(2);
%% or if x(1) position is given as an input
 %x(1) = displacementInput(1);
 %x(end) = displacementInput(2);
...
end

[Jonathan Frutschy]

@Aquatris What would you return as the state-space model dX?

[Aquatris]

You do not change dX but manually add it to your X after you solve it or you can also add to your dX with (displacementInput/dT).

Essentially what ode45 solves is

x_new = x_previous + dX*dt;

So with displacement input, you are not changing dX directly but instead you are adding something to x_previous and the resulting forces due to this displacement input appears in the dX.

[Jonathan Frutschy]

I see. Would you set F equal to zero? @Aquatris

[Jonathan Frutschy]

@Aquatris Or would you do this?

displacement_input = @(t,rd) rd.*scale.*sin(t) + epsilon.*escale; % [m]
function dX = odefun(t, rd, X, N, M, C, K,displacement_input)
    %% Definitions
    x       = X(1:N);               % position state vector
    x(1)    = x(1)+displacement_input(t,rd); 
    x(end)  = x(end)+displacement_input(t,rd);
    dx      = X(N+1:2*N);           % velocity state vector
    %% Force vector
    f_reservoir_1 = x(1);
    f_reservoir_end = x(end);
    F       = [f_reservoir_1; zeros(98,1); f_reservoir_end];
    %% Equations of Motion
    ddx     = M\(F - K*x - C*dx);
    
    %% State-space model
    dX      = [dx; ...
               ddx];
end

[Aquatris]

You can but with only this change, you will not see the displacement inputs in your X. I am not sure how you can incorporate it with ode45 solver. You might want to implement the solver your solves which is a trival for loop to mitigate the issue.

The demonstration of the issue due to using ode45:

• Time 0: x0
• Time 1: x1 = x0 + dX*dt;
• Time 2: x1_new = x1 + dispInput; -> you apply the displacement input
• Time 2: x2 = x1 + dX*dt -> in here dX is calculated with the x1_new but ode45 uses x1, not x1_new in the addition, so your x2 is actually not right