Solving 2D ordinary differential equations:

Let’s look at a problem in 2 dimensions.

Springy Pendulum

A mass on the end of a springy pendulum has a trajectory that is impossible to express analytically. The free-body diagram for the mass has forces in both the x and y directions due to 2 forces, the tension in the spring and the gravitational pull downward. The resulting acceleration of the mass is:

Where Fs = -k(L - Le), and k is the spring constant, m is the mass, L is the length of the spring, and Le is the unstretched length of the spring. Replacingand eitheror we can see that and depend only on x, y, and constants.

This time we will leave the motion in Cartesian coordinates because there is no advantage to working in polar coordinates since the radius of the motion changes as does the acceleration in all directions.

function derivs = springyPendulumODE(t, currentValues, k, m, Le)

% Springy Pendulum: (odefun)

%

% dx/dt = v_x

% dv_x/dt = -k/m(sqrt(x.^2+y.^2) – Le).*x./sqrt(x.^2+y.^2)

% dy/dt = v_y

% dv_y/dt = -g -k/m(sqrt(x.^2+y.^2) – Le).*y./sqrt(x.^2+y.^2)

%

% Our job is to set: dervis = [dx/dt; dv_x/dt; dy/dt; dv_y/dt]

% using the currentValues = [theta; omega].

g = 9.80665; % m/s^2

x = currentValues(1);

vx = currentValues(2);

y = currentValues(3);

vy = currentValues(4);

dxdt = vx;

dydt = vy;

dvxdt = -k/m*(sqrt(x.^2+y.^2)-Le).*x./sqrt(x.^2+y.^2);

dvydt = -g - k/m*(sqrt(x.^2+y.^2)-Le).*y./sqrt(x.^2+y.^2);

derivs = [dxdt; dvxdt; dydt; dvydt];


Matlab will find the solution with its ODE solver, ode45.

k = 10; % spring constant (N/m)

m = 1; % mass (kg)

Le = 0.30; % equilibrium length (m)

initialValues = [0.2, 0, -0.23, 0]; % x0, vx0, y0, vy0

timeSpan = [0 50]; % seconds

options = odeset('RelTol', 1e-10);

[t,solution] = ode45(@(t,cv) springyPendulumODE(t,cv,k,m,Le), ...

timeSpan, initialValues, options);

x = solution(:,1);

vx = solution(:,2);

y = solution(:,3);

vy = solution(:,4);

plot(x,y)

title('trajectory');

xlabel('x')

ylabel('y')

Be sure you understand these two pieces of code. There isn’t much typing, but every bit of it does very specific things.

Assignment: M10_ODE2D

1. Describe the input data and the derived data in the springy Pendulum example. Then test the code as completely as possible by setting up scenarios where you know the solution.

2. A large flying object will be affected by air resistance. The force can be approximated by with the velocity and its magnitude, and b is a constant that depends on the object and the medium.

Set up a program to calculate the acceleration, velocity, and position of a bowling ball (m = 8 kg) launched from a height of 1.5 m above the ground with an initial speed of 40 m/s at an angle of 35o. Do this over a range of values of the drag coefficient, say from b = 0 to b = 5 N-s2/m2, and plot out the trajectories (y versus x) until the ball hits the ground for each b value. (Simulate for a longer time span than is needed and then truncate the arrays that ode45 returns to y>=0.) Explain how you know the code works. (Test it.) Use find and interpolate to determine the exact position that the object hits the ground for each value of b. Make a plot of the maximum horizontal distance the ball travels as a function of b. What is a reasonable value of the drag coefficient? (You may need to Google “air resistance” to make a well reasoned estimate.)