cs249 Homework 3 Solutions

3.1)  Here is my version of interest2.m

function result = interest2 (a, p, n)
% f = interest (a, p)
% computes the minimum interest needed to reach the amount A
% after n periods, given a periodic payment of P
    i = 0.01;

    for j = 1:20
        f = func (a, p, n, i);
        fp = deriv (a, p, n, i);

        %fprintf ('%f\n', f); 
        i = i - f / fp;
        fprintf ('%f   %.8f\n', j, i);    
    end
    fprintf ('%f   %f\n', j, i); 
    result = i;
end


function f = func (a, p, n, i)
    f = p/i * (1 - (i+1)^-n) - a;
end


function fp = deriv (a, p, n, i)
    fp = -p/i^2 * (1 - (i+1)^-n) + p/i * n * (i+1)^(-n-1);
end

Notice that the escape sequence %.8f inside the printf
statements allows me to display the interest rates with
8 digits after the decimal point.

With a montly payment of $1800, we would need the interest
rate to be...

>> interest2 (640000, 1800, 360)
1.000000   -0.01959072
2.000000   -0.01644838
3.000000   -0.01323101
4.000000   -0.00994399
5.000000   -0.00664914
6.000000   -0.00356750
7.000000   -0.00121002
8.000000   -0.00011346
9.000000   0.00006500
10.000000   0.00006897
11.000000   0.00006897

The 10th and 11th interations have converged, so we assume
that the answer is close to 0.00006897, which is 0.006897%.

That is a very small interest rate, so it is probably not
feasible, in the real world, to get a mortgage with a payment
of $1800.

If we try to make a payment of $1700...

>> interest2 (640000, 1700, 360)
1.000000   -0.02198639
2.000000   -0.01889643
3.000000   -0.01573930
4.000000   -0.01250854
5.000000   -0.00922089
6.000000   -0.00597182
7.000000   -0.00305715
8.000000   -0.00105465
9.000000   -0.00032167
10.000000   -0.00024670
11.000000   -0.00024600

We converge on a negative interest rate, which means the
bank would have to pay us to take a loan.  In other words,
it is not possible.

The point of this problem is that models sometimes produce
unexpected results, like negative numbers.  In those cases,
it sometimes takes some thought to interpret the result of
the model.


3.2) The hardest part is formulating this problem as a root-finding
problem.  We are trying to minimize the amount of material in the
can, but (at first glance) it is a function of radius, height, and
the seam width.

    m = 2 pi (r + s)^2 + h 2 pi (r + s) 

But the seam width is a constant parameter, and we can write
height as a function of r (for a given volume):

    v = pi r^2 h
    h(r) = v / (pi r^2)

Plugging h(r) into m yields m as a function of a single
variable, r.  Now, we want to minimize m(r), which means finding
the roots of m'(r), the derivative.

Here is my version of can.m

function can ()
    v = 1000;
    s = 0.25;
    r = newton (v, s)
    h = v / pi / r^2
end

function result = newton (v, s)
    r = 10;

    for j = 1:10
        fprintf ('%f   %f\n', j, r); 
    
        f = func (v, s, r);
        fp = deriv (v, s, r);

        r = r - f / fp;
    end
    result = r;
end


function f = func (v, s, r)
    f = 4 * pi * (r + s) - 2 * v / r^2 - 4 * v * s / r^3;
end

function fp = deriv (v, s, r)
    fp = 4 * pi + 4 * v / r^3 + 12 * v * s / r^4;
end


The function named func evaluates m'(r).  deriv evaluates
the derivative of m'(r).

When I ran it, the result converged quickly:

>> can
1.000000   10.000000
2.000000   3.608269
3.000000   4.703636
4.000000   5.364531
5.000000   5.493282
6.000000   5.496730
7.000000   5.496733
8.000000   5.496733
9.000000   5.496733
10.000000   5.496733

r =  5.4967
h =  10.5352

This is a little squatter than I expected, but there it is.

There is one more question I wish I had asked:

If we increase the seam width, do you expect the optimal
can shape to get taller or shorter?  

Think about it and see what your instinct tells you before
you look at the answer below.

When the seam width is 0.35, the answer is 

r =  5.5248
h =  10.4283

So the can gets a little shorter and wider.