PHYS 220 — HW #4 [Solution]
The equations of motion for a golf ball in a realistic atmosphere
are given by (be careful of signs!):
- d²x/dt² = - (Fdrag_x +
Sωvy) / m
- d²y/dt² = -g - (Fdrag_y -
Sωvx) / m
The drag force is the same "quadratic" drag force we encountered
in Lab #3 and lecture. S and ω
represent a coefficient used to estimate the Magnus force on a
spinning ball and the rotation rate of that ball.
- Write out the figurative representation of the vector containing
the variables in the problem as vars = [ ... ].
Figurative:
y = [
x
y
vx
vy
]
- Using valid MATLAB code, write the single line of code for the
derivative function (if you defined an auxiliary variable for ease
of writing out that line, provide it as
well) dvarsdt = [ ... ].
MATLAB:
vx = y(3);
vy = y(4);
v = sqrt(vx^2 + vy^2);
dydt = [
vx
vy
-(b*v*vx + S*w*vy)/m
-g - (b*v*vy - S*w*vx)/m % NOTE SIGN CHANGE!
];
- List the parameters you need to pass into your derivative function.
You'll need five parameters this time:
- m → mass of golf ball
- g → acceleration due to gravity
- b → quadratic drag coefficient
- S → Magnus coefficient
- w → spin rate of ball
The following problem concerns error analysis/control. You're modelling a
non-linear oscillator that has a characteristic amplitude of
45°. You're currently running your simulation with a time step
Δto = 0.05 s. For the problem at
hand you desire a fractional error of 0.5% of the amplitude scale.
For these problems, I wrote a
short MATLAB
function to calculate the new time step. Below I show the calls
to it.
- At some point in your simulation you
estimate a local truncation error of -0.09°. Calculate
Δt1 required for the following methods:
- Euler method
new_dt(0.05,-0.09,45*0.005,2) = 0.079
- Mid-point method (rk2)
new_dt(0.05,-0.09,45*0.005,3) = 0.067
- Runge-Kutta method (rk4)
new_dt(0.05,-0.09,45*0.005,5) = 0.060
Further into the simulation you estimate a local truncation error
of 0.9°. Calculate Δt1 required for the
following methods:
- Euler method
new_dt(0.05,0.9,45*0.005,2) = 0.025
- Mid-point method (rk2)
new_dt(0.05,0.9,45*0.005,3) = 0.031
- Runge-Kutta method (rk4)
new_dt(0.05,0.9,45*0.005,5) = 0.037
Comment on the results found in 4 & 5, above.
A few things to take notice of here. In each case I
intentionally round down to a value lower than the 5-digit
MATLAB result, ensuring that I still meet the error criteria with
a smaller time step than is required. In §§4 you
see that our time step is unnecessarily small for what we need out
of the problem and that we can make it larger. Also, you can get
away with larger time steps when switching to Runge-Kutta 4 from
Euler, etc. In §§5 we have the opposite problem: too
large a time step, so we see it made smaller for each method.
Again, rk4 lets us get away with the largest time steps due to its
improved precision.