Not sure what's happening in this code and not sure where to come from here.

1 次查看（过去 30 天）
Steve Avilesmejia
Steve Avilesmejia 2020-12-2
回答： Alan Stevens 2020-12-3
clc; clear all; close all
L = 30;
m = 68.1;
c_d = 0.25;
k = 40;
gamma = 8;
g = 9.81;
h = 0.1;
t = 0: h: 50;
n = length(t);
v = zeros(2, n);
v(:,1) = [0;0];
y = t(2)
for i = 2:n
if y(1) <= L
dy = [y(2);
g - sign(y(2))*c_d*y(2)^2/m];
else
dy = [y(2);
g - sign(y(2))*c_d*y(2)^2/m - k*(y(1)-L)/m - gamma*y(2)/m];
v(:,i) = v(:, i-1) + fun(v(:,i-1), L, m, c_d, k, gamma, g)*h;
end
end

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

回答（1 个）

Alan Stevens
Alan Stevens 2020-12-3
You are not keeping track of parameters through time. Try
L = 30;
m = 68.1;
c_d = 0.25;
k = 40;
gamma = 8;
g = 9.81;
h = 0.1;
t = 0: h: 50;
n = length(t);
% You need dx/dt = v; and dv/dt = acceleration
% Simple Euler approach
x = zeros(1,numel(t)); % Allocate storage space
v = zeros(1,numel(t)); % Allocate stirage space
x(1) = 0; % Initial position
v(1) = 0; % Initial velocity
for i = 1:numel(t)-1
x(i+1) = x(i) + h*v(i); % update position
% Select acceleration
if x(i) < L
accn = g - sign(v(i))*c_d*v(i)^2/m;
else
accn = g - sign(v(i))*c_d*v(i)^2/m - k*(x(i)-L)/m - gamma*v(i)/m;
end
v(i+1) = v(i) + h*accn; % update velocity
end
% Plot position and velocity vs time
plot(t,x,t,v),grid
xlabel('time'),ylabel('Posn and Vel')
legend('Posn','Vel')

类别

在 Help Center 和 File Exchange 中查找有关 Particle & Nuclear Physics 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by