Change Input Force to Input Displacement ODE Solver
37 次查看(过去 30 天)
显示 更早的评论
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 个评论
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
回答(0 个)
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Ordinary Differential Equations 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!