COMPUTER SCIENCE ASSIGNMENT #1- MATLAB function

Assignment X, where X is the number of the corresponding assignment.
1
A MATLAB function M-file for Euler’s method (as described in class and on pages 17-19
of the textbook) for solving the differential equation
dv
dt = g −
c
m
v
(this is (1.9) on page 14) is given below. As discussed in class, this numerical method
is obtained by approximating dv
dt at time ti by v(ti+1)−v(ti)
ti+1−ti
, which results in the computed
approximation
v(tt+1) = v(ti) + h
g −
c
m
v(ti)
i
(ti+1 − ti).
To create a MATLAB function M-file, either type edit in the Command Window or
select HOME → New → Script or select HOME → New → Function (the latter gives
you a template for creating a function).
Each of these options will open a new window (an Editor window) in which to type in
the MATLAB statements for Euler’s method. Enter the following. Statements starting with
% are comments, documenting the MATLAB code.
function Euler(m,c,g,t0,v0,tn,n)
% print headings and initial conditions
fprintf(’values of t approximations v(t)\n’)
fprintf(’%8.3f’,t0),fprintf(’%19.4f\n’,v0)
% compute step size h
h=(tn-t0)/n;
% set t,v to the initial values
t=t0;
v=v0;
% compute v(t) over n time steps using Euler’s method
for i=1:n
v=v+(g-c/m*v)*h;
t=t+h;
fprintf(’%8.3f’,t),fprintf(’%19.4f\n’,v)
end
To save your M-file, select
EDITOR → Save → Save As...
At the top of this window, it should say
File name: Euler.m
Save as type: MATLAB files (*.m)
Select
Save
to save your file, and close the Editor window.
In order to use the above function, you must specify values for the 7 local parameters in
the function Euler:
2
m is the mass of the falling object
c is the drag coefficient
g is the gravity constant
t0 is the initial time, v0 is the initial velocity
tn is the final time at which the velocity is to be computed
n is the number of time steps into which [t0, tn] is divided
Thus, in the function Euler, the step size h = (tn − t0)/n is computed, and Euler’s method
is used to compute an approximation to the solution v(t) of the differential equation at the
n points (values of time)
t0 + h, t0 + 2h, t0 + 3h, t0 + 4h, . . . , t0 + nh = tn.
For example, in order to use Euler to solve the problem given in Example 1.2 on page 17
and to save your results in a file for printing, you could enter the following (in the MATLAB
Command Window):
diary filename
Euler ( 68.1 , 12.5 , 9.81, 0 , 0 , 12 , 6 )
{the desired results should appear here}
diary off
There are three questions and each subquestion (for example Q1b) is worth 2 points.
The total number of points for this assignment is 20.
Question #1. - 8 marks
(a) Create a working copy of the MATLAB function Euler in your installation of MATLAB.
DELIVERABLES: A copy of the M-FILE in your pdf.
(b) Use Euler to solve the differential equation using m = 86.2, c = 12.5 and initial
conditions v(0) = 0 on the time interval [0, 12] using 15 time steps and g = 9.81.
DELIVERABLES: Your answer should include the function call to Euler and the
resulting output.
(c) Use Euler for a falling parachutist with the same parameters as Q1b but with a gravitational constant of 3.71 (as would be the case if the parachutist was falling on Mars).
DELIVERABLES: Your answer should include the function call to Euler and the
resulting output.
3
(d) Use MATLAB to compute the relative error |εt
| in the computed approximation at
time t = 12, using the constants from Q1b. To do this, use the true (exact) solution:
v(t) = gm
c
(1 − e
− ct
m ) (1)
DELIVERABLES: A copy of the commands and output.
Note that in MATLAB e
x
is computed as exp(x) and |x| as abs(x).
Question #2 - 6 Marks.
In our mathematical model of a falling parachutist, instead of assuming that air resistance
is linearly proportional to velocity (that is, FU = −cv), you might choose to model the
upward force on the parachutist as a second-order relationship,
FU = −kv2
where k is a second-order drag coefficient. This leads to the following differential equation
dv
dt = g −
k
m
v
2
(a) Modify the MATLAB function Euler in Question 1 so that it will use Eulers method
to solve this differential equation. Use the function header
Euler2(m , k , g, t0 , v0 , tn , n)
DELIVERABLES: A copy of the M-FILE in your pdf.
(b) Use Euler2 to compute a numerical approximation to the above differential equation
using m = 73.5, k = 0.234 and initial condition v(0) = 0 on the time interval [0, 18]
using 72 time steps.
DELIVERABLES: The function call to Euler2 and the resulting output.
(c) Use the fact the exact analytic solution of this problem is
v(t) = r
gm
k
tanh r
gk
m
t
!
to compute (either in MATLAB or using your calculator) the relative error in the
computed solution at t = 18.
Question #3 - 4 Marks
The function e
−x
can be approximated by its McLaurin series expansion as follows (note
the alternating + and −):
e
−x ≈ 1 − x +
x
2
2! −
x
3
3! + . . . ±
x
n
n!
4
Alternatively, note that e
−x =
1
e
x . Thus, e
−x
can also be apporximated by 1 over the
McLaurin series expansion of e
x
. That is,
e
−x ≈
1
1 + x +
x2
2! +
x3
3! + . . . +
xn
n!
Approximate e
−2 using both approaches above for n = 1, 2, 3, 4 and 5. Note, n is the
degree of the polynomial not the number of terms. So here you use 2 terms, then 3 terms,
..., and finally 6 terms. Compare each approximation to the true value of e
−2 = 0.135335...,
using the true relative error. What conclusions can you make about the two approaches?
5