Change Input Force to Input Displacement ODE Solver

37 次查看（过去 30 天）

Jonathan Frutschy 2024-6-24，13:37

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2131446-change-input-force-to-input-displacement-ode-solver

编辑： Aquatris 2024-6-25，15:15

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?

8 个评论
显示 6更早的评论隐藏 6更早的评论

Jonathan Frutschy 2024-6-25，14:37

编辑：Jonathan Frutschy 2024-6-25，14:38

在 MATLAB Online 中打开

@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 2024-6-25，15:15

编辑：Aquatris 2024-6-25，15:15

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

请先登录，再进行评论。

请先登录，再回答此问题。