brachistochrone problem

Brachistochrone Problem PhL 10.29.16 Here I just finish the solution started by Goldstein page 36. To avoid confusion, I am going to assume the particle starts at t = 0 and y(0) = 0 and goes down to negative y (G likes that to be positive y which just makes a mess for me). We then have f = (1) y ≤ 0 Here means dy/dx. Have Maple compute the required stuff. // I did this, and I end up with this ODE = 0 (see brachisto.mws) or 2 y(x) y"(x) + y'(x)2 + 1 = 0 (2) Wolfram claims this is equivalent to [ 1 + y'(x)2] y(x) = constant (3) Proof: If I differentiate the latter, I get [ 1 + y'(x)2] y'(x) + [2y'(x) y"(x)] y(x) = 0 or [ 1 + y'(x)2] + [2 y"(x)] y(x) = 0 (4) and sure enough, this is exactly my equation (2). I have just shown that (3) true (2) true Can we go the other direction? This is what I just showed above: d { [ 1 + y'(x)2] y(x) } /dx = [ 1 + y'(x)2] + [2 y"(x)] y(x) or d { [ 1 + y'(x)2] y(x) } = { [ 1 + y'(x)2] + [2 y"(x)] y(x)}dx This is a perfect differential, so integrate this from x = x0 to some t to get [ 1 + y'(x)2] y(x) - [ 1 + y'(x0)2] y(x0) = !Syntax Error, I{ [ 1 + y'(x)2] + [2 y"(x)] y(x)}dx Now if I know that { [ 1 + y'(x)2] + [2 y"(x)] y(x)} = 0, then I know that [ 1 + y'(x)2] y(x) = [ 1 + y'(x0)2] y(x0) and this then has a specific constant. So what I know is (2) true (3) true with the specific constant shown. Comment: Be careful, if you set x0 = 0 and y(x0) = 0 you get infinite slope times 0 here. This part of the curve only applies if your final point is directly underneath your starting point. To summarize: Equation 2 2 y(x) y"(x) + y'(x)2 + 1 = 0 Equation 3 [ 1 + y'(x)2] y(x) = K Theorem 1: Eq 3 true Eq 2 true for any value of K in Eq 3 Theorem 2: Eq 2 true Eq 3 but only for K = [ 1 + y'(x0)2] y(x0) Now Eq 2 is what really comes out of the Euler-Lagrange analysis. Resume: Now I only have to solve a first order ODE [ 1 + y'(x)2] y(x) = C/2 (3) We are then supposed to try this "cycloid" solution which is the ONLY solution in fact, x = C(θ - sinθ) y = C(1 - cosθ) dx/dθ = C(1 - cosθ) dy/dθ = +Csin(θ) ≡ dy/dx = sinθ/(1-cosθ) [ 1 + y'(x)2] y(x) = [ 1 + sin2θ/(1-cosθ)2 ]C (1-cosθ) = K ?? or [(1-cosθ)2 + sin2θ ] C(1-cosθ) = K (1-cosθ)2 ?? or [(1-cosθ)2 + sin2θ ]C = K (1-cosθ) ?? or [1+cos2θ - 2cosθ + sin2θ ] C = K (1-cosθ) ?? or [2 - 2cosθ ]C = K (1-cosθ) ?? or 2C = K ?? This is true if K = C/2. So then I claim that for the above cycloid solution this equation is true, [ 1 + y'(x)2] y(x) = C/2 Now in my non-Goldstein world the y axis points up as usual, so we must have C < 0 which will make y < 0 on our trajectory. For example, we then have Remember that y axis goes up, x axis goes to the right, particle starts at origin and drops down the curve. It happens that the right side of the curve above has θ < 0 as Maple plotting shows. Where is the bottom of this curve? Look for dy/dx = 0. But above dy/dx = sinθ/(1-cosθ) so the bottom occurs at θ = -π. And at that point y = 2C. And also x = -Cπ and y = 2C. Ratio is -y/x = 2/π ≈ 2/3. Interpret 2C = the lowest value that y attains on a hump. C = (1/2) ymin < 0. Digression: The rolling wheel of radius R. Here is the rolling wheel model (it rolls to the right) Relative to the center of the wheel the marked point has coordinates x' = -Rsinθ and y' = -Rcosθ. The center of the wheel is at xc = Rθ and yc = R. Thus, the expressions for the point (x,y) are x = xc+ x' = Rθ -Rsinθ = R(θ-sinθ) y = yc + y' = R-Rcosθ = R(1-cosθ) So the dot painted on the wheel traces out the right side of the second red curve shown above, R > 0 here. Back to the problem: What then exactly is the curve? It starts at the origin, and you want it to go through some point (x,y) where x > 0 and y < 0. So consider (where we know C < 0) x = C(θ - sinθ) y = C(1-cosθ) . What do we learn from the origin point (x,y) = (0,0)? Only that θ = 0. But from the other point we x1 = C(θ - sinθ) - h1 ≡ -|y1| = C(1-cosθ) So this is two equations in two unknowns. Here hi is the drop distance > 0. Take the ratio to get x1/h1 = - (θ - sinθ)/(1-cosθ) Use this to solve for θ, then use either of the above to solve for C which we know is C < 0: So to make the brach path start at the origin and pass through the point (x,y) = (1,-2), we find that the solution angle parameter is θ = -1.40 radians (implies right side), and C = -2.41 (neg as expected). Selected Plots and Conclusions The point we want to pass through is (x,y) = (x1, -h1) where x1 and h1 are positive. Define a ratio ratio = x1/h1. If this ratio is small, you do a dive. If it is large, you dive and come back up again. Specific examples h1=1 x1=0.2 ratio = 0.2 h1=1 x1=1 ratio = 1 h1=1 x1=π/2 ratio = π/2 (zero slope at the end of the run) h1=1 x1=2 ratio = 2 (starting to come back up) h1=1 x1=4 ratio = 4 (lots of coming back up) h1=1 x1=20 ratio = 20 (coming WAY back up) h1=1 x1=100 ratio = 100 (coming WAY WAYback up) In this kind of limit, the ratio of x/2 at half way to ymax is π/2 ≈ 3/2 as noted above. Another way to say it is that x/depth = 2 (x1/2/depth) = 2 (π/2) = π. So the depth is 1/π times the distance for a long trip. I found a note on Brachistochrone trains here: https://brilliantuniverse.wordpress.com/2014/04/02/trains-and-tunnels-the-brachistochrone/ Question: What is the speed at any point on such a trip? Question: What is position at any given time? x = C(θ - sinθ) y = C(1-cosθ) . vx(t) = C(∂tθ - cosθ∂tθ) = C(θ-cosθ)(∂tθ) vy(t) = C(1-cosθ) = Csinθ(∂tθ) How do you do stuff like this when you only have a parametric form?? At some instant you know it is sliding at this angle dy/dx = sinθ/(1-cosθ) = tan(φ) We know ds = dx But this just tells you ratio of speeds. Acceleration along the curve might be this as = g sin(φ) = d2s/dt2 // could use above to get some as(θ) Then work back to velocity and then position. Day Two 10.30.16 Idea for finding θ(t): use the two equations obtained by setting ax = 0 and ay = -g by twice diffing the expressions for x and y, all in Maple. Here are the results x = C(θ - sinθ) y = C(1-cosθ) . vx = C(1 - cosθ)θ' vy = C(sinθ θ') // since C<0, θ<0,vy<0 must have θ' < 0 ax = C[(1 - cosθ)θ" + sinθ θ'2] ay = C( cosθ θ'2 + sinθ θ") Now I should be able to inject "the physics" into the problem using Newton's law: ax = 0 ay = - g // because y is up in my world, gravity pulls down STOP!!! I am using F = ma but I am neglecting the constraint force in these equations !!!! So ignore everything in blue below this point in this doc!! Resume after The two equations are then 0 = C[(1 - cosθ)θ" + sinθ θ'2] -g = C[ cosθ θ'2 + sinθ θ"] wrong The first equation tells us that (1 - cosθ)θ" + sinθ θ'2 = 0 which relates θ" to θ' which seems helpful. Try now to eliminate θ'2 : θ'2 = - [ (1 - cosθ)/sinθ ]θ" or θ" = - sinθ/(1-cosθ) * θ'2 (*) Note: for small negative θ, θ" > 0. Note: Examine ay for small negative θ: ay = C( cosθ θ'2 + sinθ θ") - + - + // cannot determine the sign of ay in this way Then the second equation says -g/C = cosθ θ'2 + sinθ θ" = - cosθ[ (1 - cosθ)/sinθ ]θ" + sinθ θ" = [ - cosθ(1 - cosθ)/sinθ + sinθ] θ" = [ - cosθ(1 - cosθ) + sin2θ] θ" / sinθ = [ - cosθ + cos2θ + sin2θ] θ" / sinθ = [ - cosθ +1] θ" / sinθ So we finally have an ODE for θ (1- cosθ) θ" + (g/C) sinθ = 0 Now use half angles 2sin2(θ/2) θ" + (g/C) 2sin(θ/2)cos(θ/2) = 0 or sin(θ/2) θ" + (g/C) cos(θ/2) = 0 or θ" + (g/C) cot(θ/2) = 0 How can I verify that I did this correctly? Here is an attempt: Last line says θ" = - θ'2 (1+cosθ)/sinθ which disagrees with my (*) above which was θ" = - sinθ/(1-cosθ) * θ'2 Is it possible that (1+cosθ)/sinθ = sinθ/(1-cosθ) ? (1+cosθ) = sin2θ/(1-cosθ) ? (1-cosθ)(1+cosθ) = sin2θ ? 1 - cos2θ = sin2θ ? yes! So we agree on that result. I then do more Maple and it seems to say ay = -C θ'2 But something is wrong here because I know C< 0 but want ay = -g, sides then have different signs. Can I get that result by hand? Back up to this point, ax = C[(1 - cosθ)θ" + sinθ θ'2] ay = C( cosθ θ'2 + sinθ θ") θ'2 = - [ (1 - cosθ)/sinθ ]θ" or θ" = - sinθ/(1-cosθ) * θ'2 (*) Use the second in ay to get ay = C( cosθ θ'2 + sinθ [ - sinθ/(1-cosθ)θ'2]) = Cθ'2( cosθ + sinθ [ - sinθ/(1-cosθ)] ) = Cθ'2( cosθ - sin2θ/(1-cosθ) ) = Cθ'2( cosθ(1-cosθ) - sin2θ )/(1-cosθ) = Cθ'2( cosθ-cos2θ - sin2θ )/(1-cosθ) = Cθ'2( cosθ-1 )/(1-cosθ) = -Cθ'2 OK, it is verified. But I still have the sign of sides problem! We now have this fairly simple ODE to deal with -Cθ'2 = ay = - g Cθ'2 = g θ'2 = (g/C) θ' = - θ(t) = - t + θ(0) Now θ(0) = left end of curve = 0, so we have θ(t) = - t as the connection between the two "parameters" . There must be some trivial way I could know this. Energy conservation says mg(-y) = (1/2)m(vx2+vy2) Above I had vx = C(1 - cosθ)θ' vy = C(sinθ θ') y = C(1-cosθ) Then energy says -mg C(1-cosθ) = (m/2)C2 [ (1-cosθ)2 + sin2θ ] θ'2 -g (1-cosθ) = (1/2)C [ (1-cosθ)2 + sin2θ ] θ'2 -g (1-cosθ) = (1/2)C [ 1 + cos2θ - 2cosθ + sin2θ ] θ'2 -g (1-cosθ) = (1/2)C [ 2 - 2cosθ ] θ'2 -g (1-cosθ) = C [ 1 - cosθ ] θ'2 -g = C θ'2 Paradox: In the equation ay = -Cθ'2 we have a sign disagreement on the two sides. 1. Why is C negative? Go back to y = C(1-cosθ). We know 1-cosθ >0 and y < 0, therefore we know that C is negative. 2. Why do I think ay = - g? The mass is being pulled down by gravity. Down is the minus y direction because that is how I have set things up. 3. Is x really positive? I have x = C(θ - sinθ) = C(-|θ| + sin|θ|) for small neg θ on the right, so x = -C(|θ| - sin|θ| ) here |θ| - sin|θ| > 0 and C < 0 to x > 0, correct. 4. Again from y = C(1-cosθ) for small θ we have y < 0 if C<0 and this is correct. Resume. The ramp has a normal force N and if we draw a little picture, we see that N + g must point along the ramp [wrong]. Total force is F = N + g and g = -g with g>0. So write F = N -g and then we have N = F + g and Nx = Fx and Ny = Fy + g Fx = Nx + gx= Nx ok Fy = Ny + gy = Ny - g ok Fy/Fx = (Ny - g)/Nx ok We then need F to point down along the curve [wrong]. Then wrong Fy/Fx = dy/dx = sinθ/(1-cosθ) // from far above - + - + - + We then have wrong (Ny - g)/Nx = sinθ/(1-cosθ) Now we can use F = ma so // F = ma = N + g Fx = max Fy = may ok ok wrong ok ay/ax = Fy/Fx = dy/dx = sinθ/(1-cosθ) ok ay = [sinθ/(1-cosθ)] ax wrong - - + + At least we now have a relationship between the two acceleration components [wrong] . Maple tells us that ax = C[(1 - cosθ)θ" + sinθ θ'2] ay = C( cosθ θ'2 + sinθ θ") // checked this with new Maple code Divide to get wrong ay/ax = [ cosθ θ'2 + sinθ θ"] / [(1 - cosθ)θ" + sinθ θ'2] = [sinθ/(1-cosθ)] Don't know if this will help, but at least it is a new fact of interest: [ cosθ θ'2 + sinθ θ"] / [(1 - cosθ)θ" + sinθ θ'2] = [sinθ/(1-cosθ)] or (1-cosθ)[ cosθ θ'2 + sinθ θ"] = [sinθ] [(1 - cosθ)θ" + sinθ θ'2] or (1-cosθ)[ cosθ θ'2 + sinθ θ"] - [sinθ] [(1 - cosθ)θ" + sinθ θ'2] = 0 or θ'2{ (1-cosθ)cosθ - sin2θ } + θ"{(1-cosθ)sinθ - [sinθ] [(1 - cosθ)] } = 0 or θ'2{ (1-cosθ)cosθ - sin2θ } = 0 or θ'2{cosθ-cos2θ - sin2θ } = 0 or θ'2{cosθ-1 } = 0 θ' = 0 which of course is a wrong result. [ and now I know why it is wrong ] Paradox of 8:15AM This says that θ'(t) = 0 which says θ(t) = constant = θ(0) = 0 so θ(t) = 0 which is total garbage. I checked the algebra in red. There is then a mental logic error somewhere in the above foot of file. Here are the possibilities of what is wrong: 1. F = ma is wrong where F is the total force on m (seems unlikely that this is wrong). 2. Fy/Fx = dy/dx is wrong. ( drew careful picture, seems right) wrong 3. ax = C[(1 - cosθ)θ" + sinθ θ'2] one or both is wrong ay = C( cosθ θ'2 + sinθ θ") (both lines copied directly from Maple and x,y correctly entered) 4. dy/dx = sinθ/(1-cosθ) is wrong (but this is verified below) 5. Algebra above is done wrong. (Maple says θ'(t) = 0 as well, see below) Look at item 4: x = C(θ - sinθ) y = C(1-cosθ) . dx = C(dθ - cosθdθ) dy = C( + sinθdθ) dy/dx = sinθ/ (1-cosθ) agrees with 4 as quoted. Let Maple do algebra of step 5. First we have this Then we do this The only way this does NOT say θ' = 0 is if the denominator vanishes. That would require -θ" + cosθ θ" - sinθ θ'2 = 0 or θ"(cosθ-1) = sinθ θ'2 θ" = [sinθ/(cosθ-1)] θ'2 But the second line above says θ"(cosθ-1) - sinθ θ'2 = 0 or θ"(1-cosθ) + sinθ θ'2 = 0 . But recall above ax = C[(1 - cosθ)θ" + sinθ θ'2] So this little option says ax= 0, but I know that is wrong because object does accelerate to the right as if slides down the ramp. How else would it pick up vx > 0 after the vertical start? Status of paradox 8:40AM. I made a list of 5 possible error sources, and checked each one and none of these sources appear to be wrong. Paradox lives on!!! Did a whiteboard session and I found the mistake. It is item 2. Although v must point along the track, in general a and therefore F do not point along the track. Think ball on string going around, or perhaps circular track with track surface pointing to the center. In that case a points to the center, it does not point along the track. It is v which points along the track. Did I make this same mistake somewhere in Lagrange doc??? Well there I had a linear track as an example in Appendix D, so in that case the force F was along the track, whew! OK, maybe this is what I want to look at : it is v which points along the track, not a. sinθ/ (1-cosθ) = dy/dx = vy/vx = C(sinθ θ') / C(1 - cosθ)θ' = sinθ/(1-cosθ) OK, this is true, but adds no new information. So here is what I know x = C(θ - sinθ) y = C(1-cosθ) . vx = C(1 - cosθ)θ' vy = C(sinθ θ') ax = C[(1 - cosθ)θ" + sinθ θ'2] ay = C( cosθ θ'2 + sinθ θ") dy/dx = vy/vx = sinθ/ (1-cosθ) How about using this extra fact from energy conservation (1/2)mv2 = mg(-y) This says v2 = -2gy or vx2 + vy2 = -2gy or C2(1 - cosθ)2θ'2 + C2(sin2θ θ'2) = -2g C(1-cosθ) or C(1 - cosθ)2θ'2 + C(sin2θ θ'2) = -2g (1-cosθ) or Cθ'2[ (1 - cosθ)2 + sin2θ] = -2g (1-cosθ) or Cθ'2[ 1 + cos2θ - 2cosθ + sin2θ] = -2g (1-cosθ) or Cθ'2[2 - 2cosθ ] = -2g (1-cosθ) or 2Cθ'2[1 - cosθ ] = -2g (1-cosθ) or Cθ'2 = -g or θ'2 = -g/C or θ' = - This seems perhaps reasonable. Then θ(t) = - t + θ(0) But θ(0) = 0, so conclude that θ(t) = - t and finally I have the connection between θ(t) and t and it is completely trivial. Marion does not mention this equation. Maple verifies this now So I can now write out everything with θ(t) = - t : x = C(θ - sinθ) y = C(1-cosθ) . vx = C(1 - cosθ)θ' vy = C(sinθ θ') ax = C[(1 - cosθ)θ" + sinθ θ'2] ay = C( cosθ θ'2 + sinθ θ") Fx = max Fy = may Nx = Fx = max Ny = Fy + g = may + g dy/dx = vy/vx = sinθ/ (1-cosθ) This is everything you ever wanted to know. But I want something to determine C from my assumed point on the trajectory which is say x1,y1 : x = C(θ - sinθ) y = C(1-cosθ) . Want to eliminate the angle. Second says y/C = 1-cosθ or cosθ = 1-y/C or θ = cos-1(1-y/C). Now for our range -2π < θ < 0 you have to be careful about the cos-1 function. But ignore for a moment and continue sinθ = - = - = - Then you get x = C( cos-1(1-y/C) - ) which you would have to solve for C in terms of x,y . Try other x2 = C2(θ - sinθ)2 y2 = C2(1 - cosθ)2 (x2+y2)/C2= (θ - sinθ)2 + (1 - cosθ)2 = θ2 + sin2θ - 2θsinθ + 1 + cos2θ - 2cosθ = θ2- 2θsinθ + 2 - 2cosθ = not very simple I think it is only numerically solvable, a "transcendental equation". Maple cannot solve it, for example, Maple can only do it numerically. Question: what is the duration of a Bracky train trip on a flat earth? Answer: θ(t) = - t so set θ = -2π to get 2π = t t = 2π But the article quoted above says Now for the special problem where y = 0, we do know that y = C(1-cosθ) = 0 cosθ = 1 θ = -2π for our situation x = C(θ - sinθ) = C(-2π -0) = 2π(-C) In this case we know that C = - (x/2π) so my result above gives t = 2π = 2π = and this agrees with the quoted result. In this case we have y = - (x/2π) (1-cosθ) and this takes its larges neg value at θ = -π so we get ymin = - (x/2π) (2) = - x/π which is about 1/3 of x H