Test 4 will cover 3 projects: Newton's Method, optimization, and integration:
using MTH229
using Roots
using Plots; plotly()
using QuadGK
Base.ctranspose(f::Function) = Roots.D(f)
tangent(f,c) = x -> f(c) + f'(c) *(x-c)
function riemann(f::Function, a::Real, b::Real, n::Int; method="right")
if method == "right"
meth = (f,l,r) -> f(r) * (r-l)
elseif method == "left"
meth = (f,l,r) -> f(l) * (r-l)
elseif method == "trapezoid"
meth = (f,l,r) -> (1/2) * (f(l) + f(r)) * (r-l)
elseif method == "simpsons"
meth = (f,l,r) -> (1/6) * (f(l) + 4*(f((l+r)/2)) + f(r)) * (r-l)
end
xs = a + (0:n) * (b-a)/n
as = [meth(f, l, r) for (l,r) in zip(xs[1:end-1], xs[2:end])]
sum(as)
end
Newton's method is an iterative method to solve $f(x) = 0$ based on the observation that if you have a good "guess" to a solution then this guess can be improved if we replace it with the solution to "the tangent line to $f$ at $x$ is $0$." Here is a picture:
f(x) = x^5 - x - 1 # x0 = 1.75 x1 = x0 - f(x0) / f'(x0) # plot(f, 1, 2, legend=false) plot!([x0, x0], [0, f(x0)]) plot!(tangent(f, x0), 1, 2) scatter!([x0, x1, newton(f, x0)], [0,0, 0])
Can you label the following on the above picture: the initial guess, the improved guess, the actual solution, the function, the tangent line?
We see the key step in Newton's method is:
$$~ x_{i+1} = x_i - f(x_i) / f'(x_i). ~$$We can compute this by hand with julia with repeated calls to x = x - f(x)/f'(x) or use one of:
newton(f, x0) where x0 is the initial guess. This function has an informative verbose=true argument to show the steps
fzero(f, x0) where x0 is the initial guess. This function is like Newton's method, but is a bit less sensitive to the quality of the initial guess.
Some examples:
Here are 3 steps of Newton's method to solve $x^5 - x - 1 = 0$ starting at $2$:
f(x) = x^5 - x - 1 x = 2 x = x - f(x)/f'(x) x = x - f(x)/f'(x) x = x - f(x)/f'(x) x, f(x)
(1.2235784352087191, 0.5190000898999054)
We see that our value of f(x) is nowhere near 0 – it is $0.519...$ – we need more steps. Here are 3 more:
x = x - f(x)/f'(x) x = x - f(x)/f'(x) x = x - f(x)/f'(x) x, f(x)
(1.1673039842671697, 4.974788514466866e-8)
Much better, our error is now $10^{-8}$, so the first 7 decimal points, $1.1673039$, are good. One more step will finish it off:
x = x - f(x)/f'(x) x, f(x)
(1.1673039782614187, 6.661338147750939e-16)
At this point, as f(x) is so close to 0, the update step of f(x) /f'(x) will basically be $0$ and the guesses will stop improving.
It took a lot of steps. We can do better if we pick more judiciosly our starting point, so we should have graphed first:
plot(f, 1, 2)
Clearly one would start closer to the answer, say at 1.2, or if we too tired to interpolate, maybe 1.5, which can be seen from the axes.
We can do this type of problem manually, as shown above, or with either of the functions:
newton(f, 2)
Or
a = fzero(f, 2) a, f(a)
(1.1673039782614187, 6.661338147750939e-16)
So, relatively easy to do, but it won't always work. We discussed three possible issues: the initial guess is too far from the answer; the first derivative at the guess is too close to $0$; or the second derivative is so big, that the tangent line approximation isn't good.
The harder problems involve solving for intersection points. For example, if $v(b) = 10*exp(-b^2)$ we can solve for a solution to $v(b) = 8$ with:
v(b) = 10 * exp(-b^2) h(b) = v(b) - 8 fzero(h, 1)
Optimization problems involve finding the largest or smallest value for some problem (largest area, least surface area, least time,...). The problems generally involve:
representing the value as a function of a single variable – usually incorporating some constraint equation
using the extreme value theorem to reduce the problem to looking either at critical points or end points
using a numeric function (e.g., fzero) to solve for the points.
The hardest part is getting started. Here is one of my least favorite optimization problems:
A farmer has 2400 ft of fencing and wants to fence off a rectangular field that borders a straight river. He needs no fence along the river. What are the dimensions of the field that has the largest area?
The two equations are: $A=x \cdot y$ and $P = 2400 = x + 2y$. The latter is the constraint – $x$ and $y$ can't be too big. We use the constraint to solve for $x$, then plug this into $A$:
A(y) = (2400 - 2y) * y # x = 2400 - 2y
A (generic function with 1 method)
The value of $y$ can range from 0 to 1200, both of which no farmer would ever consider, as these fields would have 0 area. (Why)
We plot to see:
plot(A, 0, 1200)
Here the answer is seen by symmetry, but suppose we don't see that. We are looking for the $y$ value that maximizes the graph. It is near $500$ and it is a critical point, hence solves $A'(y) = 0$:
y0 = fzero(A', 500) y0, A(y0)
(600.0, 720000.0)
Integration can be done using the fundamental theorem of calculus when an antiderivative is known. This is not always the case. If not, we can approximate the integral using the interpretation that the integral $\int_a^b f(x) dx$ has a value given by the area under the graph of $f(x)$ between $[a,b]$ when $f(x)$ is non-negative and continuous.
How to approximate the area? We discuss a few things:
the riemann function implements a Riemann sum approach or related approaches (Trapezoid and Simpson's).
the quadgk function implements a numerically more sophisticated approach (faster, more accurate)
Simple usage is straightforward:
f(x) = sin(x) riemann(f, 0, pi, 1000) # 1000 indicates how many rectangles
Or with more rectangles
riemann(f, 0, pi, 1000000)
The more rectangles, the more accuracy (at least in theory, but this will break down as soon as floating point issues arise). Other shapes give more accuracy with the same number of divisions of $[a,b]$:
riemann(f, 0, pi, 1000, method="trapezoid")
riemann(f, 0, pi, 1000, method="simpsons")
The quadgk function doesn't need to have a number of rectangles specified as it works adaptively. It also returns an estimated error:
quadgk(f, 0, pi)
(2.0, 1.7905676941154525e-12)
The project discusses two applications of the integral, let's agree to focus on just one – the arc length.
The length of the graph of $f(x)$ is given in terms of an integral related to $f$:
$$~ \int_a^b \sqrt{1 + f'(x)^2} dx. ~$$So the length of the sine curve from $0$ to $2$ is given not by quadgk(sin,0,2) – which is the area under the graph – but rather involves a related function, called dl below:
f(x) = sin(x) dl(x) = sqrt(1 + f'(x)^2) quadgk(dl, 0, 2)
(2.3516888074007873, 1.4495493694255401e-10)