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 the 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 juliabox.com 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 1: 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 2: 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 3: Newtonish methods¶

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.)

In [ ]:

Problem 4: Addition, another way¶

In [ ]:

# 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 [ ]:

x = rand(Float32, 10^7)
@time div2sum(x)
@time div2sum(x)
@time div2sum(x)
@time sum(x)
@time sum(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: