snowplow problem

[EdS2]

Quote:One day it started snowing at a heavy and steady rate. A snowplow started out at noon, going 2 miles the first hour and 1 mile the second hour. What time did it start snowing?

Eventually, I was able to solve this with high school calculus and a little algebra. But it should also be soluble by a bit of programming, perhaps by simulation, or by successive approximation, or perhaps even using SOLVE. Perhaps there's a graphical or geometric approach.

via Bob (rprosperi) via a video via a defunct blog via a 1942 textbook

The blog said (mild spoiler):

Quote:Despite being drawn from a text on the lofty subject of differential equations, this problem should be solvable by anyone who has had college-level Calculus I, or a year of high school calculus. The hardest part of the problem is not the math required to solve it, but the fact that you must make some assumption on your own about the relationship between snow and plowing speed. That is, how does the speed of a snowplow depend on the depth of the snow on the ground? The simplest reasonable assumption is this: The speed of a snowplow is inversely proportional to the depth of the snow it is plowing at each moment. Though it may be hard to believe, this assumption and the information given in the problem are enough to calculate a definite solution.

Back in college, I solved the homework problem for my mathematical modeling course by hand on paper using a little calculus, but I have often wondered how a student might approach the problem if they didn't know any calculus, but did have some ability to program. Many (if not all) “calculus problems” can be solved through fairly simple computer programming, and it turns out that this one is no exception.

(Personally, and for reasons of search, I would equally well call this the snowplough problem.)

[Albert Chan]

Let snow started at time 0, time x = noon-time, when snowplow was started.
Assumed speed of snowplow inverse-proportional to snow thickness.

First hour snowplow distance traveled = doubled of 2nd hour (due to more snow)

∫(1/t, t=x .. x+1) = 2 * ∫(1/t, t=x+1 .. x+2)

ln(x+1) - ln(x) = 2 * (ln(x+2) - ln(x+1))
3*ln(x+1) = 2*ln(x+2) + ln(x)

(x+1)³ = (x+2)² * x
x² + x - 1 = 0
(x + φ) * (x - (φ-1)) = 0

Since x > 0, x = φ-1 ≈ 0.6180 hour (37 minutes)

---

We can estimate integrals with mid-point rules.


∫(1/t, t=x .. x+1) ≈ (1/(x+1/4) + 1/(x+3/4)) / 2   // estimate with 2 rectangles
∫(1/t, t=x+1 .. x+2) ≈ 1/(x+3/2)                         // 1 rectangle OK, since the left is flatter

(x+1/2) / ((x+1/4) * (x+3/4)) = 2 * 1/(x+3/2)
(x+1/2) * (x+3/2) = 2 * (x+1/4)*(x+3/4)            // scale away denominator
x² + 2x + 3/4 = 2x² + 2x + 3/8
x² = 3/8
x = √(3/8) ≈ 0.6124 hours (37 minutes)

[EdS2]

Oh, that's nice that a rectangular approximation will do pretty well!

[Albert Chan]

(03-18-2021 08:38 PM)Albert Chan Wrote:  First hour snowplow distance traveled = doubled of 2nd hour (due to more snow)

Extending the problem to distance traveled ratio = k, k > 1

∫(1/t, t=x .. x+1) = k * ∫(1/t, t=x+1 .. x+2)

ln((x+1)/x) = k * ln((x+2)/(x+1))
ln(x/(x+1)) = -k * ln(2 - x/(x+1))

Let z = x/(x+1), 0 < z < 1

z = (2-z)^(-k)

f(z) = z - (2-z)^(-k)
f'(z) = 1 - k*(2-z)^(-k-1)

Solve by Newton's method, and we need a rough guess.
(estimate both integrals with single rectangle)

1/(x+1/2) = k /(x+3/2)
x = (3-k)/(2k-2)             -- rough guess for x
z = x/(x+1) = (3-k) / ((3-k)+(2k-2)) = (3-k)/(k+1)

Code:

function snowplow(k, verbal)
    local z = max(0,(3-k)/(k+1))
    repeat
        if verbal then print(z/(1-z)) end
        local y = (2-z)^-k
        local eps = (y-z) / (1 - k*y/(2-z))
        z = z + eps         -- newton correction
    until z == z + eps*eps
    return z/(1-z)          -- time when started snowplow
end

lua> snowplow(1.5, true)
1.4999999999999998
1.566890844473769
1.568114464087791
1.5681148594401135
1.5681148594401568       --> 94 minutes

lua> snowplow(2, true)
0.49999999999999994
0.6136363636363635
0.6180278644242209
0.6180339887380153
0.618033988749895        --> 37 minutes

[Ren]

I'm amazed at your answers, I would've thrown up my arms and said
"there's not enough information to figure that out!"

(I hated Math for many years, until too late I discovered I could've been very good with it if I had accepted it from the beginning.)

[Albert Chan]

(03-19-2021 04:46 PM)Albert Chan Wrote:  Let z = x/(x+1), 0 < z < 1

z = (2-z)^(-k)

f(z) = z - (2-z)^(-k)
f'(z) = 1 - k*(2-z)^(-k-1)

With heavy snowstorm, k is big.
Guess of z=0 is good, apply Newton's method on it:

z = 0 - f(0)/f'(0) = 2^(-k)/(1-k*2^(-k-1)) = 1/(2^k-k/2)
x = z/(1-z) = 1/(2^k-k/2-1)

With almost no snow (k≈1) , we have x → inf
Or, matching above form: x = 1/(2^k-2 - c*(k-1))

We already shown k=2 ⇒ x = φ-1 ≈ 0.618
This suggested c = 2-1/x = 2-1/(φ-1) = 2-φ ≈ 0.382

The fit is pretty good

Note: table is minutes between snow started and snowplow started, for different k's

lua> fmt = function(...) print(('%.3g\t%.3g\t%.3g'):format(...)) end
lua> for k=1.25,5,0.25 do fmt(k, 60*snowplow(k), 60/(2^k-2-0.382*(k-1))) end

1.25   212    212
1.5    94.1   94.1
1.75   55.7   55.7
2      37.1   37.1
2.25   26.3   26.3
2.5    19.5   19.5
2.75   14.8   14.8
3      11.5   11.5
3.25   9.04   9.02
3.5    7.2    7.18
3.75   5.79   5.77
4      4.69   4.67
4.25   3.82   3.8
4.5    3.12   3.11
4.75   2.57   2.56
5      2.12   2.11

[EdS2]

Great to see you having fun with this problem Albert!

Indeed Ren, it's a pity that many are put off from early encounters with mathematics. I'm a great believer in lifelong education, and it helps with this, I think. It's never too late.

[Albert Chan]

(03-19-2021 04:46 PM)Albert Chan Wrote:  Solve by Newton's method, and we need a rough guess.
(estimate both integrals with single rectangle)

1/(x+1/2) = k /(x+3/2)
x = (3-k)/(2k-2)             -- rough guess for x
z = x/(x+1) = (3-k) / ((3-k)+(2k-2)) = (3-k)/(k+1)

Playing with the formula, I discovered a simpler and better guess for x
Start with original equation:

ln((x+1)/x) = k * ln((x+2)/(x+1))
ln(1 - 1/(1+x)) = -k * ln(1 + 1/(1+x))

Let y = 1/(1+x), 0 < y < 1

ln(k) = ln(-ln(1-y)/ln(1+y)) = y + y^3/4 + 19/144*y^5 + ...

atanh(y) = y + y^3/3 + y^5/5 + ...

→ ln(k) ≈ atanh(y)                // when y is tiny (k ≈ 1)

y ≈ tanh(ln(k)) = (k-1/k) / (k+1/k) = (k²-1) / (k²+1)
x = 1/y - 1 = 2/(k²-1)       // rough guess for x

This is a better guess for x. Bonus: it has the right range, x > 0
(with this new guess x, guess z = x/(x+1) = 2/(k²+1), also very simple)

For k=1.5, new guess x=1.6 (96 minutes), old guess x=1.5 (90 minutes)
True x ≈ 1.568 (94 minutes)

For k=2.0, new guess x=2/3 (40 minutes), old guess x=1/2 (30 minutes)
True x = φ-1 ≈ 0.6180 (37 minutes)

---

Update: if we apply asinh on top of atanh, it fit ln(k) even better.

asinh(atanh(y)) = y + y^3/6 + 13/120*y^5 + ...

→ ln(k) ≈ asinh(atanh(y))

y = tanh(sinh(ln(k))) = tanh((k-1/k)/2) = (e^(k-1/k)-1)/(e^(k-1/k)+1)

x = 1/y - 1 = 2/expm1(k-1/k)

expm1(k-1/k) = 2*(k-1) + (k-1)² + (k-1)³/3 + (k-1)^4/6 - (k-1)^5/15 + ...

Keep 2 terms, x = 2/(2*(k-1)+(k-1)²) = 2/(k²-1), matching old estimate
Keep 3 terms, x = 2/(k^3/3+k-4/3) = 6/((k-1)*(k²+k+4))

This is an improvemnt over old estimate:

For k=1.5, guess x = 6/3.875 = 1.548 (93 minutes)
For k=2.0, guess x = 6/10 (36 minutes)

[Albert Chan]

(03-23-2021 05:50 PM)Albert Chan Wrote:  ln(k) = ln(-ln(1-y)/ln(1+y)) = y + y^3/4 + 19/144*y^5 + ...

Instead of getting a formula from k to x, it is better to do z = ln(k), to x.

→ y = z - z^3/4 + O(z^5) = z / (1 + z^2/4 + O(z^4)) ≈ z/(1+z^2/4)
→ x = 1/y - 1 = 1/z - 1 + z/4     // rough guess of x, from z = ln(k)

For k=1.5, guess x ≈ 1.568 (94 minutes)
For k=2.0, guess x ≈ 0.6160 (37 minutes)

This estimate matched pretty well. (Mathematica session)
Unmatched coefficients (z^5, z^7, z^9) have same sign, and similar size.

Code:

In[1]:= logk = Log[-Log[1 - y]/Log[1 + y]];

In[2]:= Series[logk, {y,0,10}]    (* z in terms of y *)

             3       5        7         9
            y    19 y    751 y    2857 y        11
Out[2]= y + -- + ----- + ------ + ------- + O[y]
            4     144     8640     44800

In[3]:= InverseSeries[%] /. y->z  (* y in terms of z *)

             3    5       7         9
            z    z    91 z    3379 z        11
Out[3]= z - -- + -- - ----- + ------- + O[z]
            4    18   8640    1814400

In[4]:= Series[z/(1+z^2/4), {z,0,10}] (* Pade[2,2] *)

             3    5    7    9
            z    z    z    z         11
Out[4]= z - -- + -- - -- + --- + O[z]
            4    16   64   256

Update: more terms, x in terms in z = ln(k)

x = 1/z - 1 + z/4 + z^3/144 - 7/4320*z^5 - 73/3628800*z^7 + 11/1451520*z^9 + ...