HW3 (1)
pdf
keyboard_arrow_up
School
University of Florida *
*We aren’t endorsed by this school
Course
3344
Subject
English
Date
Apr 3, 2024
Type
Pages
4
Uploaded by MinisterMooseMaster790
EGM 3344
Summer 2023
Introduction to Numerical Methods of Engineering Analysis
University of Florida
Mechanical and Aerospace Engineering
Homework 3: These Are the Roots You’re Looking For
Due: Saturday, June 10, 2023, 11:59pm.
Instructions:
1. Submit your solutions as a PDF document on Canvas.
2.
Include any and all plots in the PDF.
3. Also submit your individual .m files.
4. Do not submit a zipped file.
5. You may
not
use symbolic variables.
6. You
may
use function handles.
7. Remember to use a superfluous amount of comments!
Problem 1:
Following the discussion in problem 5.10 (oxygen concentration) from the textbook,
write a MATLAB script called
hw3_p1.m
to use the bisection method to determine the
temperature in degrees Celsius for which
σ
sf
= 14 mg/L. Use 0
◦
C and 32
◦
C as the
initial guesses and approximate absolute error as the stopping criterion with
ϵ
s
= 0
.
5,
where
ϵ
s
is called the “tolerance.” The full stopping criterion is: stop if
error < ϵ
s
.
For your convenience, the equation for the saturation concentration of dissolved oxygen
in freshwater is:
log
σ
sf
=
−
139
.
34411 +
1
.
575701
×
10
5
T
a
−
6
.
642308
×
10
7
T
2
a
+
1
.
243800
×
10
10
T
3
a
−
8
.
621949
×
10
11
T
4
a
,
where
σ
sf
is the saturation concentration of dissolved oxygen in freshwater at 1 atm
(mg/L) and
T
a
is the absolute temperature (K). Note the units
1
!
Discuss your results
—how many iterations of the bisection method were required?
Could you have predicted the number of iterations
a priori
2
?
Problem 2:
Write a MATLAB script called
hw3_p2.m
to find the root of
f
(
x
) = 0
.
1
e
x
x
−
1 using
the bisection method. Solve the problem six times with the following combinations of
initial brackets and stopping criteria:
Initial bracket:
1
This is not just some pretentious professorial “units are important—you have to label units.” The units affect your
code.
2
This is Latin for “from before.”
It means “before getting more information.”
In this context, it means “before
running the algorithm.”
Jonathan Brooks
1
EGM 3344
Summer 2023
i) Use
x
ℓ
= 0 and
x
u
= 2.
ii) Use
x
ℓ
=
−
3 and
x
u
= 4.
Stopping criterion:
i) Use
ϵ
a
≤
10%.
ii) Use
ϵ
a
≤
1%.
iii) Use
ϵ
a
≤
10
−
6
%.
Use for loop
s
to solve the problem for each initial-bracket/stopping-criterion pair, and
have MATLAB display your answer in the Command Window (with the initial bracket,
stopping criterion threshold, number of iterations, and root). Also include these results
in your PDF.
Discuss your results:
Do you notice any trends?
Problem 3:
Write a MATLAB script called
hw3_p3.m
to find the root of
f
(
x
) = 0
.
1
e
x
x
−
1 using the
false-position method. Solve the problem six times with the following combinations of
initial brackets and stopping criteria:
Initial bracket:
i) Use
x
ℓ
= 0 and
x
u
= 2.
ii) Use
x
ℓ
=
−
3 and
x
u
= 4.
Stopping criterion:
i) Use
ϵ
a
≤
10%.
ii) Use
ϵ
a
≤
1%.
iii) Use
ϵ
a
≤
10
−
6
%.
Use for loop
s
to solve the problem for each initial-bracket/stopping-criterion pair, and
have MATLAB display your answer in the Command Window (with the initial bracket,
stopping criterion threshold, number of iterations, and root). Also include these results
in your PDF.
Discuss your results:
Is there a reason that some of the scenarios required signifi-
cantly more iterations?
Hint:
It might be useful to plot the function on the initial interval to get a feel for
how the method might work in these scenarios.
Problem 4:
Write a MATLAB script called
hw3_p4.m
to find the root of
f
(
x
) = 0
.
1
e
x
x
−
1 using
the Newton-Raphson method for three different initial guesses: i)
x
0
= 1, ii)
x
0
=
−
1,
and iii)
x
0
=
−
2.
Jonathan Brooks
2
EGM 3344
Summer 2023
Use a for loop to solve the problem for each initial guess, and have MATLAB dis-
play your answer with the initial guess and the number of iterations in the Command
Window
3
. Also include these results in your PDF.
Discuss your results:
Why did the Newton-Raphson method work sometimes and
fail other times?
Hint:
Plot the function to gain intuition.
Problem 5:
Write a MATLAB script called
hw3_p5.m
to find the root of
f
(
x
) = 0
.
1
e
x
x
−
1 using the
secant method with two sets of initial guesses: i) (
x
0
, x
1
) = (
−
1
,
1) and ii) (
x
0
, x
1
) =
(1
,
−
1).
Use a for loop to solve the problem for each set of initial guesses, and have MAT-
LAB display your answer with the initial guesses and the number of iterations in the
Command Window
4
. Also include these results in your PDF.
Discuss your results:
Why does the algorithm either fail or take a very long time
for one set of initial guesses?
Hint:
Use debug mode and watch the
x
and
f
(
x
) values. One will be very unlike the
others.
Extra hint:
Think about your discussion for problem 3.
Problem 6:
With our rocket successfully launched, we need to determine where its target—Mars—
is. This problem was solved by Johannes Kepler by supposing that Mars has a circular
orbit instead of an elliptical one and then correcting for the actual orbit’s eccentricity.
Kepler derived the following equation (which is appropriately named Kepler’s equation):
2
πt
T
=
E
(
t
)
−
e
·
sin(
E
(
t
))
,
where
t
is time,
T
is the orbital period,
E
(
t
) is an angle (in radians) called the eccentric
anomaly, and
e
is the eccentricity of the orbit. The eccentric anomaly is the angular
position of the planet when projected onto the supposed circular orbit and changes
with time. Unfortunately, this equation cannot be solved for
E
(
t
) analytically, so we
must rely on numerical methods.
We would like the rocket to arrive at Mars in 360 days. Write a set of MATLAB scripts
called
hw3_p6a.m
,
hw3_p6b.m
,
hw3_p6c.m
, and
hw3_p6d.m
to determine the eccentric
anomaly for
t
= 360 days four times:
a) with the bisection method;
b) with the false-position method;
c) with the Newton-Raphson method;
3
Note that I am not giving you a stopping criterion. You, as the engineer, get to choose it such that you feel confident
and dignified in your answer. A dignified engineer does not use computation time or a preset number of iterations as
the stopping criterion.
4
Note that I am not giving you a stopping criterion. You, as the engineer, get to choose it such that you feel confident
and dignified in your answer. A dignified engineer does not use computation time or a preset number of iterations as
the stopping criterion.
Jonathan Brooks
3
Your preview ends here
Eager to read complete document? Join bartleby learn and gain access to the full version
- Access to all documents
- Unlimited textbook solutions
- 24/7 expert homework help
EGM 3344
Summer 2023
d) with the secant method.
Use
T
= 687 days for Mars’ orbital period and
e
= 0
.
0934 for the orbit’s eccentricity.
Use approximate relative error as the stopping criterion, with
ϵ
s
= 10
−
6
% for each
part
5
. Display your answers in the Command Window, and include them in your PDF,
as well.
Discuss your results:
Were any algorithms clearly better or worse than the others?
5 bonus points:
Write a MATLAB script called
hw3_p6e.m
to solve the problem
using fixed-point iteration.
Bonus Problem 7:
10 points
As mentioned in problem 6, Kepler would project a satellite’s position onto a circular orbit and
determine the eccentric anomaly—the angle formed by the satellite’s projected position, the center
of the orbit, and perigee. Once the eccentric anomaly is known, trigonometry can be used to obtain
formulas for the satellite’s
x
- and
y
-coordinates, where the positive
x
direction is the direction of
perigee and the positive
y
direction is the direction corresponding to
E
(
t
) =
π/
2, both measured
from the focus of the orbit:
x
=
a
(cos(
E
(
t
))
−
e
)
,
y
=
b
sin(
E
(
t
))
,
where
a
and
b
are the semi-major and semi-minor axes of the orbit, respectively.
Choose your favorite root-finding method of the five in problem 6, and write a MATLAB script called
hw3_p7.m
to use it to solve for the eccentric anomaly for every
6
value of
t
through one full orbit.
Then calculate and plot Mars’ orbit (
y
vs.
x
). Use 1
.
524 AU for Mars’ orbit’s semi-major axis along
with the eccentricity and orbital period from problem 6.
Also plot
E
(
t
) vs.
t
in a separate plot. Remember axis labels!
5
Note that I am not giving you initial brackets or guesses. You choose these as the engineer, but give each method a
reasonable chance to find the root, i.e., do not give one method a really good guess and another an absolutely terrible
guess.
6
Numerically speaking, that is.
Jonathan Brooks
4