Predicting distribution of random variable Y=AX+B from the knowledge about distribution of random variable X: a case study of projectile motion
Let us first simulate projectile motion without any
inaccuracy.
close all
clear
figure
g=9.81; %gravity
V=5; %m/s
projection speed
theta=pi/3; %projection angle
x(1)=0;%initial
x
y(1)=0;%initial
y
k=2;
for t=0.001:0.001:1
x(k)=V*cos(theta)*t+x(1);
y(k)=-0.5*g*t^2+V*sin(theta)*t+y(1);
k=k+1;
end
plot(x,y,'.')
xlim([0 3])
ylim([-1 2])
hold on
quiver(x(1),y(1),V*cos(theta),V*sin(theta),'AutoScaleFactor',0.2)
xlabel('x')
ylabel('y')
plot(0,0,'o','MarkerSize',20,'MarkerFaceColor','r')
Now we simulate the case that the initalization has
error:
The projectile is projected with velocity of V at an angle
from horizon. The projectile
should be projected from origion but the mechanism that places the projectile at the origion is
not accurate. Its error
in x direction follows a
Gaussian distribution with mean
and standard deviation of
(sigma_x) and its error
in y direction follows a
Gaussian distribution with mean
and standard deviation of
(sigma_x). These errors
are independent from one another. The following simulation code allows you to investigate the
statistical distribution of the location of the projectile after 1 second from projection.
figure
M=500;%number of
trials
g=9.81; %gravity
V=5; %m/s
projection speed
theta=pi/3; %projection angle
sigma_x=0.15; %standard deviation of the noise in x direction
sigma_y=0.05; %standard deviation of the noise in y direction
To simulate the noisy starts we use the following
commands
- x(i,1)=sigma_x*randn;%initial x
- y(i,1)=sigma_y*randn;%initial y
for i=1:M
x(i,1)=sigma_x*randn;%initial x
y(i,1)=sigma_y*randn;%initial y
k=2;
for t=0.001:0.001:1
x(i,k)=V*cos(theta)*t+x(i,1);
y(i,k)=-0.5*g*t^2+V*sin(theta)*t+y(i,1);
k=k+1;
end
subplot(1,3,1)
plot(x(:,1),y(:,1),'o')
title('Initial
conditions')
xlim([-4*sigma_x
4*sigma_x])
ylim([-4*sigma_y
4*sigma_y])
axis square
xlabel('x(0)')
ylabel('y(0)')
subplot(1,3,2)
title('Projectile
trajectory at each trial')
plot(x(i,:),y(i,:),'.')
hold on
end
xlim([0 3])
ylim([-1 2])
axis square
xlabel('x')
ylabel('y')
subplot(1,3,3)
plot(x(:,end),y(:,end),'bo')
title('Final
position after 1 sec')
hold on
axis square
xlim([2 3])
ylim([-0.9 -0.2])
You can predice the distribusion of the location of the
projectile after 1 seconds from launch using the the equations of motion and how the initial
placement erorr enters these equations.
figure
mu = [V*cos(theta) -0.5*g+V*sin(theta)];
% found from analyzing equations of motion
Sigma = diag([sigma_x^2,sigma_y^2]);
% found from analyzing equations of motion
rng('default') % For
reproducibility
R = mvnrnd(mu,Sigma,M);
plot(R(:,1),R(:,2),'r+') %these are points that
we predict from our analysis
hold on
axis square
plot(x(:,end),y(:,end),'bo') %these are points that
the simulaiton of 500 trial cases give us
axis square
xlim([2 3])
ylim([-0.9 -0.2])
xlabel('x(1)')
ylabel('y(1)')
We can see that the points predicted from the
distirbution model we obtained from our analysis closely capture the out come of 500 trials we
conducted, attesting that our statistical model is accurate.
Next, we simulate the case that the projection speed is
faulty
figure
clear x y
M=500;%number of
trials
g=9.81; %gravity
V_command=5; %m/s
projection speed
theta=pi/3; %projection angle
sigma_v=0.2;
e=sigma_v*randn(1000,1);
for i=1:M
x(i,1)=0;%initial x
y(i,1)=0;%initial y
V=V_command+e(i);
k=2;
for t=0.001:0.001:1
x(i,k)=V*cos(theta)*t+x(i,1);
y(i,k)=-0.5*g*t^2+V*sin(theta)*t+y(i,1);
k=k+1;
end
subplot(1,2,1)
plot(x(i,:),y(i,:),'.')
title('Projectile
trajectory at each trial')
hold on
end
xlim([0 3])
ylim([-1 2])
axis square
hold on
xlabel('x')
ylabel('y')
subplot(1,2,2)
plot(x(:,end),y(:,end),'o','MarkerSize',8)
axis square
xlim([2 3])
ylim([-1.1 -0.1])
hold on
You can predice the distribusion of the location of the
projectile after 1 seconds from launch using the the equations of motion and how the initial
placement erorr enters these equations.
figure
mu = [V_command*cos(theta)
-0.5*9.81+V_command*sin(theta)];
Sigma = [cos(theta)
sin(theta)]'*[cos(theta) sin(theta)]*sigma_v^2;
rng('default') % For
reproducibility
R = mvnrnd(mu,Sigma,M);
plot(R(:,1),R(:,2),'r+')
hold on
plot(x(:,end),y(:,end),'o','MarkerSize',8)
axis square
xlim([2 3])
ylim([-1.1 -0.1])
@Solmaz Kia, UCI, May 10th, 2023
-----
Solmaz Kia (she/her), Ph.D.,
Associate Professor,
UCI Inclusive Excellence Professor,
Mechanical and Aerospace Eng. Dept.,
Computer Science Dept. (joint appointment),
University of California Irvine,
solmaz.eng.uci.edu
email: solmaz@uci.edu