MATLAB: Law of motion through numerical aproximation

law of motionnumerical aproximationode

A ball of mass m=0.12 kg is a t rest at the origin of the X axis.
At the moment t0=0 a force parallel to the X axis begins to act upon the ball.
The force varies in time by the expression:
F(t)=6*sin(4*pi*t-pi/6)-0.5.
Numerically aproximate the law of motion x(t) of the ball and plot it for the time interval t=[0,10].
m=0.12 
initial_time=0;
final_time=10;
N=1001;
t=linspace(initial_time,final_time,N);
dt=t(2)-t(1);
x(1)=0;
x(2)=0;
for i=2:N-1
    x(i+1)=2*x(i)-x(i-1)+dt*dt*(6*sin(4*pi*t(i)-pi/6)-0.5);
end;
plot(t,x);
m=0.12 ;
initial_time=0;
final_time=10;
N=1001;
t=linspace(initial_time,final_time,N);
dt=t(2)-t(1);
The ball at t0 is at position 0 on the X axis, so x0=0
there is no initial speed so x1=x0=0
x(1)=0;
x(2)=0;
for i=2:N-1
    x(i+1)=2*x(i)-x(i-1)+dt*dt*(6*sin(4*pi*t(i)-pi/6)-0.5)/m;
end
plot(t,x);
I do not think I got the correct solution to this problem, it seems a little bit weird that the ball is only moving in one direction, shouldn't it oscilate araound an center point given that the function is dependant of sin?

Best Answer

Here's a basic integration of the system, along with plots:
% mass, initial time, final time
m = 0.12;
t0 = 0;
tf = 3;
% number of time intervals and resulting spacing dt
N = 1001;
t = linspace(t0,tf,N);
dt = t(2) - t(1);
% vectors for storage of position, velocity, and acceleration
x = zeros(1,N);
xdot = zeros(1,N);
xdoubledot = zeros(1,N);
%initial conditions for position, velocity, and acceleration
x(1) = 0;
xdot(1) = 0;
xdoubledot(1) = (6*sin(4*pi*t(1)-pi/6)-0.5)/m;
% numerical integration
for i = 2:1:N
    x(i) = x(i-1) + xdot(i-1)*dt;
    xdot(i) = xdot(i-1) + xdoubledot(i-1)*dt;
    xdoubledot(i) = (6*sin(4*pi*t(i)-pi/6)-0.5)/m;
end
figure(1)
plot(t,xdoubledot)
title('acceleration')
figure(2)
plot(t,xdot)
title('velocity')
figure(3)
plot(t,x)
title('position')
I left the final time at 3 seconds so you can get a better picture of what's going on. I encourage you to play around with the values of t0 and tf a little bit.
Notice that you have a constant term (-0.5) added to a sinusoidal term. What do you expect the results of this "offset" to be?
Related Question