Give your solution to Trefethen problem 13.2(c) below. Note that you can use the Insert menu (or ctrl-m b) to insert new code cells as needed.
If we are solving the quadratic equation $x^2 - 2bx + c = 0$, everyone knows the quadratic formula $x_\pm = b \pm \sqrt{b^2 - c}$. However, consider the accuracy of the $x_-$ root of this equation as $c$ gets smaller:
b = 1
c = logspace(-20,-1,50) # 50 values from 10^-20 to 10^-1, equally spaced on a log scale
x = b - sqrt(b^2 - c) # x computed in the default double precision (~15 digits)
set_bigfloat_precision(256) # 256 bits, about 80 decimal digits
xbig = b - sqrt(b^2 - big(c)) # x computed in "bigfloat" precision (~80 digits)
err = float64(abs(x - xbig) ./ abs(xbig)) # relative error, converted to double-precision to plot
using PyPlot
loglog(c, err, "bo-")
xlabel(L"c")
ylabel(L"relative error in $x$")
title(L"standard quadratic formula $x = b - \sqrt{b^2 - c}$")
Now, here is a more accurate way to compute $x_-$:
Suppose we are trying to find a root $f(x)=0$ and we are given ways to compute $f$, $f'$, and $f''$. You will design an algorithm to take advantage of this, and apply it to $f(x)=x^3-1$.
(You should be able to base your code on the Newton square-root example from class.)
# Sum x[first:last]. This function works, but is a little slower than we would like.
function div2sum(x, first=1, last=length(x))
n = last - first + 1;
if n < 2
s = zero(eltype(x))
for i = first:last
s += x[i]
end
return s
else
mid = div(first + last, 2) # find middle as (first+last)/2, rounding down
return div2sum(x, first, mid) + div2sum(x, mid+1, last)
end
end
# check its accuracy for a set logarithmically spaced n's. Since div2sum is slow,
# we won't go to very large n or use too many points
N = round(Int, logspace(1,7,50)) # 50 points from 10¹ to 10⁷
err = Float64[]
for n in N
x = rand(Float32, n)
xdouble = float64(x)
push!(err, abs(div2sum(x) - sum(xdouble)) / abs(sum(xdouble)))
end
loglog(N, err, "bo-")
title("simple div2sum")
xlabel("number of summands")
ylabel("relative error")
grid()
Time it vs. the built-in sum
(which is also written in Julia):
x = rand(Float32, 10^7)
@time div2sum(x)
@time sum(x)
You should notice that it's pretty darn slow compared to sum
, although in an absolute sense it is pretty good. Make it faster: