18.335 Problem Set 1

This notebook accompanies the first problem set posted on the 18.335 web page, and is here to get you started with your own Julia computations.

Download this notebook (a pset1.ipynb file) by right-clicking the download link at the upper-right to Save As a file, and then drag this file into your Jupyter dashboard to upload it (e.g. on Binder or in a local installation).

Modify it as needed, then choose Print Preview from the "File" menu and print to a PDF file to submit electronically.

Problem 2: Floating point

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.

In [ ]:

Problem 3: Funny functions

part (a)

Compute $(|x|^4 + |y|^4)^{1/4}$:

In [ ]:
L4(x,y) = ...

Some tests:

In [ ]:
L4(1e-100,0.0)
In [ ]:
L4(1e+100,0.0)

part (b)

Compute $\cot(x) - \cot(x + y)$:

In [ ]:
cotdiff(x,y) = ...

Some tests:

In [ ]:
cotdiff(1.0, 1e-20)
In [ ]:
cotdiff(big(1.0), big(1e-20)) # compute in BigFloat precision

Problem 4: Addition, another way

In [7]:
# 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, 10 .^ range(1,7,length=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

using PyPlot
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):

In [6]:
x = rand(Float32, 10^7)
@time div2sum(x)
@time div2sum(x)
@time div2sum(x)
@time sum(x)
@time sum(x)
@time sum(x)
  0.047769 seconds (1 allocation: 16 bytes)
  0.046905 seconds (1 allocation: 16 bytes)
  0.044345 seconds (1 allocation: 16 bytes)
  0.002188 seconds (1 allocation: 16 bytes)
  0.002432 seconds (1 allocation: 16 bytes)
  0.002175 seconds (1 allocation: 16 bytes)
4.997705f6

You should notice that it's pretty darn slow compared to sum, although in an absolute sense it is pretty good. Make it faster:

In [ ]:
function fast_div2sum(x, first=1, last=length(x))
    # ???
end