How to plot by using result data?
6 次查看(过去 30 天)
显示 更早的评论
Dear Sir or Madam,
I am new to matlab programming I would like to know how to plot position_x , position_y , velocity x , velocity y by applying result data. Please show me how to perform the plot as I would like to see them.
Hopefully and Thank you very much for your help.
n_part=10;
kn=5;
kt=2/7*kn;
m=0.3;
g=9.81;
timestep=50;
v_init=0.5;
% initial of position, velocity, angular velocity
position_x(:,1:n_part)=rand(size(1:n_part));
position_y(:,1:n_part)=rand(size(1:n_part));
theta=2*pi*rand(size(position_x));
velocity_x=v_init*sin(theta);
velocity_y=v_init*cos(theta);
w_x=zeros(size(position_x));
w_y=zeros(size(position_y));
w_half_x=zeros(size(velocity_x));
w_half_y=zeros(size(velocity_y));
w_x_dot=zeros(size(velocity_x));
w_y_dot=zeros(size(velocity_y));
v_half=rand(size(velocity_x));
v_half_x=rand(size(velocity_x));
v_half_y=rand(size(velocity_y));
dt=0.001;
rad(1:n_part)=0.5;
acc_x=zeros(size(position_x));
acc_y=zeros(size(position_y));
acc_i=[0,0];
acc_j=[0,0];
w_i=[0,0];
w_j=[0,0];
Fn=[0,0];Fn_x=[0,0];Fn_y=[0,0];
for t=1:1:4
for i=1:n_part
v_half_x(t+1,i)=velocity_x(t,i)+0.5*dt*acc_x(t,i);
v_half_y(t+1,i)=velocity_y(t,i)+0.5*dt*acc_y(t,i);
position_x(t+1,i)=position_x(t,i)+v_half_x(t+1,i)*dt;
position_y(t+1,i)=position_y(t,i)+v_half_y(t+1,i)*dt;
% update of w_ part
w_half_x(t+1,i)=w_x(t,i)+dt*0.5*dt*w_x_dot(t,i);
w_half_y(t+1,i)=w_y(t,i)+dt*0.5*dt*w_y_dot(t,i);
j=1;
for j=i+1
if j<=n_part
lx=position_x(i)-position_x(j);
ly=position_y(i)-position_y(j);
vx=velocity_x(i)-velocity_x(j);
vy=velocity_y(i)-velocity_y(j);
root_xy=sqrt(lx.^2+ly.^2);
R_i=rad(i);
R_j=rad(j);
delta=(R_i+R_j)-root_xy;
if delta>0
N=[lx;ly];
n_vector=N/norm(N);
% relative vector to find tangent vector;
vij1=[vx;vy];
w_i=[w_x(i);w_y(i)];
w_j=[w_x(j);w_y(j)];
vij2=vij1+(R_i.*w_i+R_j.*w_j).*n_vector;
v_i_j=-n_vector.*(n_vector.*vij2);
% tangent vector
tan=v_i_j/norm(v_i_j);
Fn=kn.*delta;
Ft=kt.*v_i_j;
F_total=Fn.*n_vector+Ft.*tan+m*g;
Fn_x=Fn_x-F_total(1,:);
Fn_y=Fn_y+F_total(2,:);
%Momemtum of particles
M_i=R_i.*n_vector.*Ft;
M_j=R_j.*n_vector.*Ft;
% same radius for movement of inertia of circle
I=pi*R_i/4;
w_i=0;
w_j=0;
w_i=w_i+M_i'/I;
w_j=w_i+M_j'/I;
w_x_dot(t,i)=w_x_dot(t,i)+w_i(1);
w_x_dot(t,i)=w_x_dot(t,i)+w_j(1);
w_y_dot(t,i)=w_y_dot(t,i)+w_i(2);
w_y_dot(t,i)=w_y_dot(t,i)+w_j(2);
w_x_dot(t+1,:)=w_x_dot(t,:);
w_y_dot(t+1,:)=w_y_dot(t,:);
acc_i=0;
acc_i=acc_i+Fn_x./m;
acc_j=0;
acc_j=acc_j+Fn_y./m;
acc_x(t,i)=acc_x(t,i)+acc_i(1);
acc_y(t,i)=acc_y(t,i)+acc_i(2);
acc_x(t,j)=acc_x(t,j)+acc_j(1);
acc_y(t,j)=acc_y(t,j)+acc_j(2);
acc_x(t+1,:)=acc_x(t,:);
acc_y(t+1,:)=acc_y(t,:);
end
end
end
velocity_x(t+1,i)=v_half_x(t+1,i)+0.5*dt*acc_x(t+1,i);
velocity_y(t+1,i)=v_half_y(t+1,i)+0.5*dt*acc_y(t+1,i);
w_x(t+1,i)=w_half_x(t+1,i)+0.5*dt*w_x_dot(t+1,i);
w_y(t+1,i)=w_half_y(t+1,i)+0.5*dt*w_y_dot(t+1,i);
end
end
% my intension is using rad, position_x, position_y to create circle as
% particles and using velocity_ x and velocity _y to move them.
% position_x (t,i); position_y (t,i); velocity_x (t,i); velocity_y (t,i);
% rad=0.5;
0 个评论
回答(1 个)
另请参阅
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!