Using ode45 to solve 2nd order nonlinear ODE

Hello, I'm trying to figure out how to use the ode45 solver on MATLAB to solve for the problem below.

Given a ball with the mass of M = 0.25kg is launched up into the air from the ground at a initial velocity of v0 = 50m/s at an angle of 30 degrees to the surface of the ground. The force acting on the ball is F = Cv^(3/2) and the differential equations are x'' = -(C/M)*x'*sqrt[(x')^2)+(y')^2)] and y'' = (-g)-(C/M)*y'*sqrt[(x')^2)+(y')^2)]. The values of g = 9.81m/s^2 and C = 0.03kg/(ms)^(1/2). Find the time of flight and range of the ball.

I've broken down the 2 2nd order ODE's into 4 1st order ODE

U1 = x

U2 = x'

U3 = y

U4 = y'

U1' = U2

U2' = -(C/M)*U2*sqrt[(U2)^2+(U4)^2]

U3' = U4

U4' = (-g)-(C/M)*U4*sqrt[(U2)^2+(U4)^2]

How do I use ode45 solver to solve for problem? I need the plot for the data vs time along with the phase space (x vs x').

[Sam Chak]

Hi @Fue Xiong

Can you share the link and examples on the ode45 solver so that I can read them and help to solve your problem?

[Fue Xiong]

Sorry for the confusion. MATLAB has a function called ode45() that solves differential equations. Here is the link, https://www.mathworks.com/help/matlab/ref/ode45.html#bu00_4l_sep_shared-y0

What I got so far

U1 = x (this is displacement)

U2 = x' (this is velocity)

U3 = y (this is displacement)

U4 = y' (this is velocity)

U1' = U2

U2' = -(C/M)*U2*sqrt[(U2)^2+(U4)^2]

U3' = U4

U4' = (-g)-(C/M)*U4*sqrt[(U2)^2+(U4)^2]

This is my code so far

[t x]=ode45(@FinalEP1,[0 10],[0 50 0 50]);

%Plot results

figure(1)

plot(t,x(:,1),t,x(:,2),t,x(:,3),t,x(:,4))

legend('x1','x2','x3','x4')

xlabel('Time t')

ylabel('Displacement and Velocity')

-----------------------------------

Function definition

function [x_primes] = FinalEP1(t,x)

m = 0.25;

C = .03;

g = 9.81;

x_primes = [x(2); (-C/m)*x(2)*sqrt((x(2)^2)+(x(4)^2)); x(4); -g-(C/m)*x(4)*sqrt((x(2)^2)+(x(4)^2))]

Results

I am so far able to plot the data but I am unsure how to interpret the results. I believe x1 and x3 are the displacement for the ball. x2 and x4 is the velocity with initial of 50m/s. The question asks me to find the time of flight and range of the ball.

Answer 1

[Sam Chak]

Hi @Fue Xiong

You actually need to formulate the simultaneous problem first

function F = simulprob(x)

F(1) = sqrt(x(1)^2 + x(2)^2) - 50;

F(2) = atan(x(2)/x(1)) - 30*pi/180;

and then solve it

fun = @simulprob;
x0 = [40, 20];
x = fsolve(fun, x0)

to obtain initial velocities

and

With the information, now you can numerically solve the ODE properly.

tspan = 0:1e-4:3;               % simulation time
init = [0 25*sqrt(3) 0 25];     % initial values
[t, y] = ode45(@FinalEP1, tspan, init);

and plot it to see how the states evolve

figure(1)
plot(t, y, 'linewidth', 1.5)
grid on
xlabel('Time, t [sec]')
title('Time responses of the System')
legend({'$x(t)$', '$v_{x}$', '$y(t)$', '$v_{y}$'}, 'Interpreter', 'latex', 'location', 'best')

Since the Ball is assumed to be fired from a flat ground, so it does not make sense if the vertical range in y has negative values because . So, you need to find the time t when y crosses zero, and this is how you can do it:

yA = sign(y(:,3));
yB = diff(yA);
v = find(yB)
time_of_flight = t(v(2))
X = y(:,1);
Y = y(:,3);
range_of_ball = X(v(2))

Now you can plot the Trajectory of the Ball:

figure(2)
plot(X, Y, 'linewidth', 1.5)
xlim([0 range_of_ball])
grid on
xlabel('horizontal range in x')
ylabel('vertical range in y')
title('Trajectory of the Ball')