Exercise 10 - Using MATLAB to solve an ODE

Practice

Try making the radioactive decay script work with the addition of plotting the true answer for comparison to the calculated answer. (The link to the example program is here) Make a number of runs with different step size to see the improvement of the calculation as the step size is decreased.

Calculation

We will develop a script to solve a first order ordinary differential equation by the "Euler" method. This exercise leads you through the necessary steps.

The equation we will solve is that for the velocity of a falling ball including air resistance:

(NOTE: equation has been corrected effective 10/27/09; the correction has no impact on the answer since m=1)

The exact solution is:

(NOTE: equation has been corrected effective 10/27/09; the correction has no impact on the answer since m=1)

where V₀ is the initial velocity. For this problem we will assume our ball is at rest initially, so V₀=0. The gravitational acceleration g is the traditional 9.8 m/s/s. We will use a drag coefficient, b, of 0.1 N/(m/s). The mass of our object will be assumed to be 1 kg.

We will write a MATLAB script to solve this using the iterative method described in class and display the results.

1. Write a script which calculates the exact solution of v(t) for an array of times, t (see above).

• Have the script span a range of time from 0 to 100. Make the user input the step size for the t array. To get started, try a step size of 2
• calculate the exact solution array, v, which corresponds in size to the t array. Our initial velocity is 0 (at rest); g = 9.8 m/s/s; b = 0.1 N/(m/s); m=1 kg.
• plot the answer

2. Now create another array to hold the computed solution for the set of times, t, and do the computation of the solution array using the Euler method.

• You will have to initialize the first element of the computed solution array to the initial value of V₀, corresponding to our initial condition of the ball at rest.
• Then, each time through the for loop, you will be computing the solution at the next time step defined by the t array. Think about when this will need to be ended in your loop.
• Plot the results of your calculation on the same graph as the calculation of the exact solution for comparison.

3. If you used the suggested values for time step, then the solutions above do not agree. Try using smaller and bigger time steps to see the effect of improving the approximation. When you have achieved a good solution, save the plot as a JPEG file and submit both the plot and the script used to make the plot to the directory: ~schloerb/ph281/username

Extra Credit

Write a script in which the air resistance is proportional to velocity squared rather than just linear in velocity. For a specific case where the air resistance force is equal for these two laws at a velocity of 10 m/s (i.e. bv = cv² when v=10 so if b=0.1 above then c=0.01 in the new air resistance force formula) compare the behavior of the falling ball under the two laws.

Go to Ph281 Home

Go to Ph281 Topics