|
Parameterizing Motion along a Curve
by Shawn Hedman, Florida Southern College, shedman@flsouthern.edu
Introduction
We animate a point mass moving along a parabola, a helix, and a piecewise defined combination of these two curves.
The animations model the motion of a bead sliding along a wire under a couple of assumptions.
| • |
We make the terrestrial assumption that a mass of m kg weighs 9.8m N.
|
| • |
We make the unrealistic assumption that the motion of the bead is not retarded by friction, air resistence, or other forces. That is, we assume the total energy of the system remains constant.
|
Initialization
Part 1 A point mass on a parabola.
|
|
(1) |
Suppose a point mass of m kg is released from the point (0,9) with an initial velocity of 0 m/s.
As the bead falls along the curve, it loses potential energy and gains kinetic energy.
| • |
Kinetic energy is given by: KE =  where  m/s is the speed of the point mass.
|
The motion of the point mass from time t=a to t=b minimizes ![(∫)[a]^(b)(KE-PE)dt](/view.aspx?SI=130465/430192/05992adfccbe2d2d48cb281e677fe310.gif) .
The difference  is called the Lagrangian and is denoted by L:
|
|
(2) |
To minimize this integral, the point mass will move fastest at points where the kinetic energy is greatest and slowest at points where the potential energy is greatest.
Because we are assuming the total energy remains constant:  for some C.
Initially, we have:  , and so  .
In general, when the point mass is located at  we have:
|
|
(3) |
Solving for  yields:
|
|
(4) |
Our goal is to parameterize the curve as ![r1(t) = [x1(t), y1(t)]](/view.aspx?SI=130465/430192/210ad52599e7e447543466a164dfb729.gif) with the required speed.
The velocity of this parameterization is ![v(t) = [diff(x1(t), t), diff(y1(t), t)]](/view.aspx?SI=130465/430192/541283bbed4b7fa50b910aaddaac82df.gif) .
The speed v is the magnitude of the velocity  
|
|
(5) |
Squaring both sides: 
|
|
(6) |
Note: If we set (4) equal to (6) we get a nonlinear first order ordinary differential equation. Because  this differential equation can be expressed in terms of  (or in terms of  ). The initial value problem with initial condition  (or  ) has constant solution  (or  ). Because v has been canceled out of the equation, we cannot give the object a nonzero initial velocity. This is not helpful.
To find the required parameterization ![r1(t) = [x1(t), y1(t)]](/view.aspx?SI=130465/430192/1be65d54be936c97383f1f4bc68450e0.gif) , we use the fact that the functions  and  minimize ![(∫)[a]^(b)L dt](/view.aspx?SI=130465/430192/fab8f13f74be0bc6d7f7f06ab892263d.gif) where L is the Lagrangian:
|
|
(7) |
To minimize L, we first express L in terms of a single variable.
|
|
(8) |
|
|
(9) |
|
|
(10) |
 Substituting these quantitites into (7) gives L in terms of  :
|
|
(11) |
To minimize L, the above expression must satisfy the Euler-Lagrange equation:
d/dt (∂L/∂X)=∂L/∂x.
In this equation, X represents dx/dt, but is treated as an independent variable. To compute the left and right sides of this equation, we express the Lagrangian as:
|
|
(12) |
Because m occurs on both sides of the Euler-Lagrange equation, it cancels out. For this reason, we have omitted it from the definition of L1. We have also replaced g with its numerical value.
|
|
(13) |
|
|
(14) |
EL1R is the right side of the Euler-Lagrange equation. We call the left side EL1L.
|
|
(15) |
EL1L is the derivative of the previous output with respect to t.
|
|
(16) |
|
|
(17) |
The Euler-Lagrange equation (EL) is a second order nonlinear ordinary differential equation for x1. We solve this with the initial conditions given by x=0 and v=0:
|
 else _xout := evalf(x_rkf45) end if; _dat := Array(1..4, {(1) = proc (_xin) local _xout, _dtbl, _dat, _vmap, _x0, _y0, _val, _dig, _n, _ne, _nd, _nv, _pars, _ini, _par, _i, _j, _k, _src; option `Copyright (c) 2002 by Waterloo Maple Inc. All rights reserved.`; table( [( "complex" ) = false ] ) _xout := _xin; _pars := []; _dtbl := array( 1 .. 4, [( 1 ) = (array( 1 .. 18, [( 1 ) = (datatype = float[8], order = C_order, storage = rectangular), ( 2 ) = (datatype = float[8], order = C_order, storage = rectangular), ( 3 ) = ([0, 0, 0, Array(1..0, 1..2, {}, datatype = float[8], order = C_order)]), ( 5 ) = (Array(1..20, {(1) = .0, (2) = 0.10e-5, (3) = .0, (4) = 0.500001e-14, (5) = .0, (6) = 0.31762478565669594e-2, (7) = .0, (8) = 0.10e-5, (9) = .0, (10) = .0, (11) = .0, (12) = .0, (13) = 1.0, (14) = .0, (15) = .5, (16) = .0, (17) = 1.0, (18) = 1.0, (19) = .0, (20) = .0}, datatype = float[8], order = C_order)), ( 4 ) = (Array(1..34, {(1) = 2, (2) = 2, (3) = 0, (4) = 0, (5) = 0, (6) = 0, (7) = 1, (8) = 0, (9) = 0, (10) = 0, (11) = 0, (12) = 0, (13) = 0, (14) = 0, (15) = 0, (16) = 0, (17) = 0, (18) = 1, (19) = 30000, (20) = 0, (21) = 0, (22) = 1, (23) = 4, (24) = 0, (25) = 1, (26) = 15, (27) = 1, (28) = 1, (29) = 1, (30) = 3, (31) = 3, (32) = 0, (33) = 1, (34) = 0}, datatype = integer[4])), ( 7 ) = ([Array(1..4, 1..7, {(1, 1) = .0, (1, 2) = .203125, (1, 3) = .3046875, (1, 4) = .75, (1, 5) = .8125, (1, 6) = .40625, (1, 7) = .8125, (2, 1) = 0.6378173828125e-1, (2, 2) = .0, (2, 3) = .279296875, (2, 4) = .27237892150878906, (2, 5) = -0.9686851501464844e-1, (2, 6) = 0.1956939697265625e-1, (2, 7) = .5381584167480469, (3, 1) = 0.31890869140625e-1, (3, 2) = .0, (3, 3) = -.34375, (3, 4) = -.335235595703125, (3, 5) = .2296142578125, (3, 6) = .41748046875, (3, 7) = 11.480712890625, (4, 1) = 0.9710520505905151e-1, (4, 2) = .0, (4, 3) = .40350341796875, (4, 4) = 0.20297467708587646e-1, (4, 5) = -0.6054282188415527e-2, (4, 6) = -0.4770040512084961e-1, (4, 7) = .77858567237854}, datatype = float[8], order = C_order), Array(1..6, 1..6, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (1, 5) = .0, (1, 6) = 1.0, (2, 1) = .25, (2, 2) = .0, (2, 3) = .0, (2, 4) = .0, (2, 5) = .0, (2, 6) = 1.0, (3, 1) = .1875, (3, 2) = .5625, (3, 3) = .0, (3, 4) = .0, (3, 5) = .0, (3, 6) = 2.0, (4, 1) = .23583984375, (4, 2) = -.87890625, (4, 3) = .890625, (4, 4) = .0, (4, 5) = .0, (4, 6) = .2681884765625, (5, 1) = .1272735595703125, (5, 2) = -.5009765625, (5, 3) = .44921875, (5, 4) = -0.128936767578125e-1, (5, 5) = .0, (5, 6) = 0.626220703125e-1, (6, 1) = -0.927734375e-1, (6, 2) = .626220703125, (6, 3) = -.4326171875, (6, 4) = .1418304443359375, (6, 5) = -0.861053466796875e-1, (6, 6) = .3131103515625}, datatype = float[8], order = C_order), Array(1..6, {(1) = .0, (2) = .386, (3) = .21, (4) = .63, (5) = 1.0, (6) = 1.0}, datatype = float[8], order = C_order), Array(1..6, {(1) = .25, (2) = -.1043, (3) = .1035, (4) = -0.362e-1, (5) = .0, (6) = .0}, datatype = float[8], order = C_order), Array(1..6, 1..5, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (1, 5) = .0, (2, 1) = 1.544, (2, 2) = .0, (2, 3) = .0, (2, 4) = .0, (2, 5) = .0, (3, 1) = .9466785280815533, (3, 2) = .25570116989825814, (3, 3) = .0, (3, 4) = .0, (3, 5) = .0, (4, 1) = 3.3148251870684886, (4, 2) = 2.896124015972123, (4, 3) = .9986419139977808, (4, 4) = .0, (4, 5) = .0, (5, 1) = 1.2212245092262748, (5, 2) = 6.019134481287752, (5, 3) = 12.537083329320874, (5, 4) = -.687886036105895, (5, 5) = .0, (6, 1) = 1.2212245092262748, (6, 2) = 6.019134481287752, (6, 3) = 12.537083329320874, (6, 4) = -.687886036105895, (6, 5) = 1.0}, datatype = float[8], order = C_order), Array(1..6, 1..5, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (1, 5) = .0, (2, 1) = -5.6688, (2, 2) = .0, (2, 3) = .0, (2, 4) = .0, (2, 5) = .0, (3, 1) = -2.4300933568337584, (3, 2) = -.20635991570891224, (3, 3) = .0, (3, 4) = .0, (3, 5) = .0, (4, 1) = -.10735290581452621, (4, 2) = -9.594562251021896, (4, 3) = -20.470286148096154, (4, 4) = .0, (4, 5) = .0, (5, 1) = 7.496443313968615, (5, 2) = -10.246804314641219, (5, 3) = -33.99990352819906, (5, 4) = 11.708908932061595, (5, 5) = .0, (6, 1) = 8.083246795922411, (6, 2) = -7.981132988062785, (6, 3) = -31.52159432874373, (6, 4) = 16.319305431231363, (6, 5) = -6.0588182388340535}, datatype = float[8], order = C_order), Array(1..3, 1..5, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (1, 5) = .0, (2, 1) = 10.126235083446911, (2, 2) = -7.487995877607633, (2, 3) = -34.800918615557414, (2, 4) = -7.9927717075687275, (2, 5) = 1.0251377232956207, (3, 1) = -.6762803392806898, (3, 2) = 6.087714651678606, (3, 3) = 16.43084320892463, (3, 4) = 24.767225114183653, (3, 5) = -6.5943891257167815}, datatype = float[8], order = C_order)]), ( 6 ) = (Array(1..2, {(1) = .0, (2) = .0}, datatype = float[8], order = C_order)), ( 10 ) = ([proc (N, X, Y, YP) option `[Y[1] = x1(t), Y[2] = diff(x1(t),t)]`; YP[2] := (-4*(Y[1]-3)*Y[2]^2-19.6*Y[1]+58.8)/(1+4*(Y[1]-3)^2); YP[1] := Y[2]; 0 end proc, -1, 0, 0, 0, 0]), ( 11 ) = (Array(1..6, 0..2, {(1, 1) = .0, (1, 2) = .0, (2, 0) = .0, (2, 1) = .0, (2, 2) = .0, (3, 0) = .0, (3, 1) = .0, (3, 2) = .0, (4, 0) = .0, (4, 1) = .0, (4, 2) = .0, (5, 0) = .0, (5, 1) = .0, (5, 2) = .0, (6, 0) = .0, (6, 1) = .0, (6, 2) = .0}, datatype = float[8], order = C_order)), ( 8 ) = ([Array(1..2, {(1) = .0, (2) = .0}, datatype = float[8], order = C_order), Array(1..2, {(1) = .0, (2) = .0}, datatype = float[8], order = C_order), Array(1..2, {(1) = .0, (2) = 1.5891891891891892}, datatype = float[8], order = C_order)]), ( 9 ) = ([Array(1..2, {(1) = .1, (2) = .1}, datatype = float[8], order = C_order), Array(1..2, {(1) = .0, (2) = .0}, datatype = float[8], order = C_order), Array(1..2, {(1) = .0, (2) = .0}, datatype = float[8], order = C_order), Array(1..2, {(1) = .0, (2) = .0}, datatype = float[8], order = C_order), Array(1..2, {(1) = .0, (2) = .0}, datatype = float[8], order = C_order), Array(1..2, 1..2, {(1, 1) = .0, (1, 2) = .0, (2, 1) = .0, (2, 2) = .0}, datatype = float[8], order = C_order), Array(1..2, 1..2, {(1, 1) = .0, (1, 2) = .0, (2, 1) = .0, (2, 2) = .0}, datatype = float[8], order = C_order), Array(1..2, 1..6, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (1, 5) = .0, (1, 6) = .0, (2, 1) = .0, (2, 2) = .0, (2, 3) = .0, (2, 4) = .0, (2, 5) = .0, (2, 6) = .0}, datatype = float[8], order = C_order), Array(1..2, {(1) = 0, (2) = 0}, datatype = integer[4]), Array(1..2, {(1) = .0, (2) = .0}, datatype = float[8], order = C_order), Array(1..2, {(1) = .0, (2) = .0}, datatype = float[8], order = C_order), Array(1..2, {(1) = .0, (2) = .0}, datatype = float[8], order = C_order), Array(1..2, {(1) = .0, (2) = .0}, datatype = float[8], order = C_order)]), ( 15 ) = ("rkf45"), ( 14 ) = ([0, 0]), ( 13 ) = (), ( 12 ) = (), ( 16 ) = ([0, 0, 0, 0, []]), ( 17 ) = ([proc (N, X, Y, YP) option `[Y[1] = x1(t), Y[2] = diff(x1(t),t)]`; YP[2] := (-4*(Y[1]-3)*Y[2]^2-19.6*Y[1]+58.8)/(1+4*(Y[1]-3)^2); YP[1] := Y[2]; 0 end proc, -1, 0, 0, 0, 0]), ( 18 ) = ([]) ] )) ] ); _y0 := Array(0..2, {(1) = 0., (2) = 0.}); _vmap := array( 1 .. 2, [( 1 ) = (1), ( 2 ) = (2) ] ); _x0 := _dtbl[1][5][5]; _n := _dtbl[1][4][1]; _ne := _dtbl[1][4][3]; _nd := _dtbl[1][4][4]; _nv := _dtbl[1][4][16]; if not type(_xout, 'numeric') then if member(_xout, ["start", "left", "right"]) then if _Env_smart_dsolve_numeric = true or _dtbl[1][4][10] = 1 then if _xout = "left" then if type(_dtbl[2], 'table') then return _dtbl[2][5][1] end if elif _xout = "right" then if type(_dtbl[3], 'table') then return _dtbl[3][5][1] end if end if end if; return _dtbl[1][5][5] elif _xout = "method" then return _dtbl[1][15] elif _xout = "storage" then return evalb(_dtbl[1][4][10] = 1) elif _xout = "leftdata" then if not type(_dtbl[2], 'array') then return NULL else return eval(_dtbl[2]) end if elif _xout = "rightdata" then if not type(_dtbl[3], 'array') then return NULL else return eval(_dtbl[3]) end if elif _xout = "enginedata" then return eval(_dtbl[1]) elif _xout = "enginereset" then _dtbl[2] := evaln(_dtbl[2]); _dtbl[3] := evaln(_dtbl[3]); return NULL elif _xout = "initial" then return procname(_y0[0]) elif _xout = "laxtol" then return _dtbl[`if`(member(_dtbl[4], {2, 3}), _dtbl[4], 1)][5][18] elif _xout = "numfun" then return `if`(member(_dtbl[4], {2, 3}), _dtbl[_dtbl[4]][4][18], 0) elif _xout = "parameters" then return [seq(_y0[_n+_i], _i = 1 .. nops(_pars))] elif _xout = "initial_and_parameters" then return procname(_y0[0]), [seq(_y0[_n+_i], _i = 1 .. nops(_pars))] elif _xout = "last" then if _dtbl[4] <> 2 and _dtbl[4] <> 3 or _x0-_dtbl[_dtbl[4]][5][1] = 0. then error "no information is available on last computed point" else _xout := _dtbl[_dtbl[4]][5][1] end if elif _xout = "function" then if _dtbl[1][4][33]-2. = 0 then return eval(_dtbl[1][10], 1) else return eval(_dtbl[1][10][1], 1) end if elif _xout = "map" then return copy(_vmap) elif type(_xin, `=`) and type(rhs(_xin), 'list') and member(lhs(_xin), {"initial", "parameters", "initial_and_parameters"}) then _ini, _par := [], []; if lhs(_xin) = "initial" then _ini := rhs(_xin) elif lhs(_xin) = "parameters" then _par := rhs(_xin) elif select(type, rhs(_xin), `=`) <> [] then _par, _ini := selectremove(type, rhs(_xin), `=`) elif nops(rhs(_xin)) < nops(_pars)+1 then error "insufficient data for specification of initial and parameters" else _par := rhs(_xin)[-nops(_pars) .. -1]; _ini := rhs(_xin)[1 .. -nops(_pars)-1] end if; _xout := lhs(_xout); if _par <> [] then `dsolve/numeric/process_parameters`(_n, _pars, _par, _y0) end if; if _ini <> [] then `dsolve/numeric/process_initial`(_n-_ne, _ini, _y0, _pars, _vmap) end if; `dsolve/numeric/SC/reinitialize`(_dtbl, _y0, _n, procname, _pars); if _Env_smart_dsolve_numeric = true and type(_y0[0], 'numeric') and _dtbl[1][4][10] <> 1 then procname("right") := _y0[0]; procname("left") := _y0[0] end if; if _xout = "initial" then return [_y0[0], seq(_y0[_vmap[_i]], _i = 1 .. _n-_ne)] elif _xout = "parameters" then return [seq(_y0[_n+_i], _i = 1 .. nops(_pars))] else return [_y0[0], seq(_y0[_vmap[_i]], _i = 1 .. _n-_ne)], [seq(_y0[_n+_i], _i = 1 .. nops(_pars))] end if elif _xin = "eventstop" then if _nv = 0 then error "this solution has no events" end if; _i := _dtbl[4]; if _i <> 2 and _i <> 3 then return 0 end if; if _dtbl[_i][4][10] = 1 and assigned(_dtbl[5-_i]) and _dtbl[_i][4][9] < 10 and 10 <= _dtbl[5-_i][4][9] then _i := 5-_i; _dtbl[4] := _i; _j := round(_dtbl[_i][4][17]); return round(_dtbl[_i][3][1][_j, 1]) elif 10 <= _dtbl[_i][4][9] then _j := round(_dtbl[_i][4][17]); return round(_dtbl[_i][3][1][_j, 1]) else return 0 end if elif _xin = "eventstatus" then if _nv = 0 then error "this solution has no events" end if; _i := [selectremove(proc (a) options operator, arrow; _dtbl[1][3][1][a, 7] = 1 end proc, {seq(_j, _j = 1 .. round(_dtbl[1][3][1][_nv+1, 1]))})]; return ':-enabled' = _i[1], ':-disabled' = _i[2] elif _xin = "eventclear" then if _nv = 0 then error "this solution has no events" end if; _i := _dtbl[4]; if _i <> 2 and _i <> 3 then error "no events to clear" end if; if _dtbl[_i][4][10] = 1 and assigned(_dtbl[5-_i]) and _dtbl[_i][4][9] < 10 and 10 < _dtbl[5-_i][4][9] then _dtbl[4] := 5-_i; _i := 5-_i end if; if _dtbl[_i][4][9] < 10 then error "no events to clear" else _j := _dtbl[_i][4][9]-10; if irem(round(_dtbl[_i][3][1][_j, 4]), 2) = 1 then error "retriggerable events cannot be cleared" end if; _j := round(_dtbl[_i][3][1][_j, 1]); for _k to _nv do if _dtbl[_i][3][1][_k, 1] = _j then if _dtbl[_i][3][1][_k, 2] = 3 then error "range events cannot be cleared" end if; _dtbl[_i][3][1][_k, 8] := _dtbl[_i][3][1][_nv+1, 8] end if end do; _dtbl[_i][4][17] := 0; _dtbl[_i][4][9] := 0; if _dtbl[1][4][10] = 1 then if _i = 2 then try procname(procname("left")) catch: end try else try procname(procname("right")) catch: end try end if end if end if; return elif type(_xin, `=`) and member(lhs(_xin), {"eventenable", "eventdisable"}) then if _nv = 0 then error "this solution has no events" end if; if type(rhs(_xin), {('set')('posint'), ('list')('posint')}) then _i := {op(rhs(_xin))} elif type(rhs(_xin), 'posint') then _i := {rhs(_xin)} else error "event identifiers must be integers in the range 1..%1", round(_dtbl[1][3][1][_nv+1, 1]) end if; if select(proc (a) options operator, arrow; _nv < a end proc, _i) <> {} then error "event identifiers must be integers in the range 1..%1", round(_dtbl[1][3][1][_nv+1, 1]) end if; _k := {}; for _j to _nv do if member(round(_dtbl[1][3][1][_j, 1]), _i) then _k := `union`(_k, {_j}) end if end do; _i := _k; if lhs(_xin) = "eventdisable" then _dtbl[4] := 0; _j := [evalb(assigned(_dtbl[2]) and member(_dtbl[2][4][17], _i)), evalb(assigned(_dtbl[3]) and member(_dtbl[3][4][17], _i))]; for _k in _i do _dtbl[1][3][1][_k, 7] := 0; if assigned(_dtbl[2]) then _dtbl[2][3][1][_k, 7] := 0 end if; if assigned(_dtbl[3]) then _dtbl[3][3][1][_k, 7] := 0 end if end do; if _j[1] then for _k to _nv+1 do if _k <= _nv and not type(_dtbl[2][3][4][_k, 1], 'undefined') then userinfo(3, {'events', 'eventreset'}, `reinit #2, event code `, _k, ` to defined init `, _dtbl[2][3][4][_k, 1]); _dtbl[2][3][1][_k, 8] := _dtbl[2][3][4][_k, 1] elif _dtbl[2][3][1][_k, 2] = 0 and irem(iquo(round(_dtbl[2][3][1][_k, 4]), 2), 2) = 0 then userinfo(3, {'events', 'eventreset'}, `reinit #2, event code `, _k, ` to initial init `, _x0); _dtbl[2][3][1][_k, 8] := _x0 else userinfo(3, {'events', 'eventreset'}, `reinit #2, event code `, _k, ` to fireinitial init `, _x0-1); _dtbl[2][3][1][_k, 8] := _x0-1 end if end do; _dtbl[2][4][17] := 0; _dtbl[2][4][9] := 0; if _dtbl[1][4][10] = 1 then procname(procname("left")) end if end if; if _j[2] then for _k to _nv+1 do if _k <= _nv and not type(_dtbl[3][3][4][_k, 2], 'undefined') then userinfo(3, {'events', 'eventreset'}, `reinit #3, event code `, _k, ` to defined init `, _dtbl[3][3][4][_k, 2]); _dtbl[3][3][1][_k, 8] := _dtbl[3][3][4][_k, 2] elif _dtbl[3][3][1][_k, 2] = 0 and irem(iquo(round(_dtbl[3][3][1][_k, 4]), 2), 2) = 0 then userinfo(3, {'events', 'eventreset'}, `reinit #3, event code `, _k, ` to initial init `, _x0); _dtbl[3][3][1][_k, 8] := _x0 else userinfo(3, {'events', 'eventreset'}, `reinit #3, event code `, _k, ` to fireinitial init `, _x0+1); _dtbl[3][3][1][_k, 8] := _x0+1 end if end do; _dtbl[3][4][17] := 0; _dtbl[3][4][9] := 0; if _dtbl[1][4][10] = 1 then procname(procname("right")) end if end if else for _k in _i do _dtbl[1][3][1][_k, 7] := 1 end do; _dtbl[2] := evaln(_dtbl[2]); _dtbl[3] := evaln(_dtbl[3]); _dtbl[4] := 0; if _dtbl[1][4][10] = 1 then if _x0 <= procname("right") then try procname(procname("right")) catch: end try end if; if procname("left") <= _x0 then try procname(procname("left")) catch: end try end if end if end if; return elif type(_xin, `=`) and lhs(_xin) = "eventfired" then if not type(rhs(_xin), 'list') then error "'eventfired' must be specified as a list" end if; if _nv = 0 then error "this solution has no events" end if; if _dtbl[4] <> 2 and _dtbl[4] <> 3 then error "'direction' must be set prior to calling/setting 'eventfired'" end if; _i := _dtbl[4]; _val := NULL; if not assigned(_EnvEventRetriggerWarned) then _EnvEventRetriggerWarned := false end if; for _k in rhs(_xin) do if type(_k, 'integer') then _src := _k elif type(_k, 'integer' = 'anything') and type(evalf(rhs(_k)), 'numeric') then _k := lhs(_k) = evalf[max(Digits, 18)](rhs(_k)); _src := lhs(_k) else error "'eventfired' entry is not valid: %1", _k end if; if _src < 1 or round(_dtbl[1][3][1][_nv+1, 1]) < _src then error "event identifiers must be integers in the range 1..%1", round(_dtbl[1][3][1][_nv+1, 1]) end if; _src := {seq(`if`(_dtbl[1][3][1][_j, 1]-_src = 0., _j, NULL), _j = 1 .. _nv)}; if nops(_src) <> 1 then error "'eventfired' can only be set/queried for root-finding events and time/interval events" end if; _src := _src[1]; if _dtbl[1][3][1][_src, 2] <> 0. and _dtbl[1][3][1][_src, 2]-2. <> 0. then error "'eventfired' can only be set/queried for root-finding events and time/interval events" elif irem(round(_dtbl[1][3][1][_src, 4]), 2) = 1 then if _EnvEventRetriggerWarned = false then WARNING(`'eventfired' has no effect on events that retrigger`) end if; _EnvEventRetriggerWarned := true end if; if type(_dtbl[_i][3][4][_src, _i-1], 'undefined') or _i = 2 and _dtbl[2][3][1][_src, 8] < _dtbl[2][3][4][_src, 1] or _i = 3 and _dtbl[3][3][4][_src, 2] < _dtbl[3][3][1][_src, 8] then _val := _val, _dtbl[_i][3][1][_src, 8] else _val := _val, _dtbl[_i][3][4][_src, _i-1] end if; if type(_k, `=`) then userinfo(3, {'events', 'eventreset'}, `manual set event code `, _src, ` to value `, rhs(_k)); _dtbl[_i][3][1][_src, 8] := rhs(_k); _dtbl[_i][3][4][_src, _i-1] := rhs(_k) end if end do; return [_val] elif type(_xin, `=`) and lhs(_xin) = "direction" then if not member(rhs(_xin), {-1, 1, ':-left', ':-right'}) then error "'direction' must be specified as either '1' or 'right' (positive) or '-1' or 'left' (negative)" end if; _src := `if`(_dtbl[4] = 2, -1, `if`(_dtbl[4] = 3, 1, undefined)); _i := `if`(member(rhs(_xin), {1, ':-right'}), 3, 2); _dtbl[4] := _i; _dtbl[_i] := `dsolve/numeric/SC/IVPdcopy`(_dtbl[1], `if`(assigned(_dtbl[_i]), _dtbl[_i], NULL)); if 0 < _nv then for _j to _nv+1 do if _j <= _nv and not type(_dtbl[_i][3][4][_j, _i-1], 'undefined') then userinfo(3, {'events', 'eventreset'}, `reinit #4, event code `, _j, ` to defined init `, _dtbl[_i][3][4][_j, _i-1]); _dtbl[_i][3][1][_j, 8] := _dtbl[_i][3][4][_j, _i-1] elif _dtbl[_i][3][1][_j, 2] = 0 and irem(iquo(round(_dtbl[_i][3][1][_j, 4]), 2), 2) = 0 then userinfo(3, {'events', 'eventreset'}, `reinit #4, event code `, _j, ` to initial init `, _x0); _dtbl[_i][3][1][_j, 8] := _x0 else userinfo(3, {'events', 'eventreset'}, `reinit #4, event code `, _j, ` to fireinitial init `, _x0-2*_i+5.0); _dtbl[_i][3][1][_j, 8] := _x0-2*_i+5.0 end if end do end if; return _src else return "procname" end if end if; if _xout = _x0 then return [_x0, seq(evalf(_dtbl[1][6][_vmap[_i]]), _i = 1 .. _n-_ne)] end if; _i := `if`(_x0 <= _xout, 3, 2); if _xin = "last" and 0 < _dtbl[_i][4][9] and _dtbl[_i][4][9] < 10 then _dat := eval(_dtbl[_i], 2); _j := _dat[4][20]; return [_dat[11][_j, 0], seq(_dat[11][_j, _vmap[_i]], _i = 1 .. _n-_ne-_nd), seq(_dat[8][1][_vmap[_i]], _i = _n-_ne-_nd+1 .. _n-_ne)] end if; if not type(_dtbl[_i], 'array') then _dtbl[_i] := `dsolve/numeric/SC/IVPdcopy`(_dtbl[1], `if`(assigned(_dtbl[_i]), _dtbl[_i], NULL)); if 0 < _nv then for _j to _nv+1 do if _j <= _nv and not type(_dtbl[_i][3][4][_j, _i-1], 'undefined') then userinfo(3, {'events', 'eventreset'}, `reinit #5, event code `, _j, ` to defined init `, _dtbl[_i][3][4][_j, _i-1]); _dtbl[_i][3][1][_j, 8] := _dtbl[_i][3][4][_j, _i-1] elif _dtbl[_i][3][1][_j, 2] = 0 and irem(iquo(round(_dtbl[_i][3][1][_j, 4]), 2), 2) = 0 then userinfo(3, {'events', 'eventreset'}, `reinit #5, event code `, _j, ` to initial init `, _x0); _dtbl[_i][3][1][_j, 8] := _x0 else userinfo(3, {'events', 'eventreset'}, `reinit #5, event code `, _j, ` to fireinitial init `, _x0-2*_i+5.0); _dtbl[_i][3][1][_j, 8] := _x0-2*_i+5.0 end if end do end if end if; if _xin <> "last" then if 0 < 0 then if `dsolve/numeric/checkglobals`(op(_dtbl[1][14]), _pars, _n, _y0) then `dsolve/numeric/SC/reinitialize`(_dtbl, _y0, _n, procname, _pars, _i) end if end if; if _dtbl[1][4][7] = 0 then error "parameters must be initialized before solution can be computed" end if end if; _dat := eval(_dtbl[_i], 2); _dtbl[4] := _i; try _src := `dsolve/numeric/SC/IVPrun`(_dat, _xout) catch: userinfo(2, `dsolve/debug`, print(`Exception in solnproc:`, [lastexception][2 .. -1])); error end try; if _src = 0 and 10 < _dat[4][9] then _val := _dat[3][1][_nv+1, 8] else _val := _dat[11][_dat[4][20], 0] end if; if _src <> 0 or _dat[4][9] <= 0 then _dtbl[1][5][1] := _xout else _dtbl[1][5][1] := _val end if; if _i = 3 and _val < _xout then Rounding := -infinity; if _dat[4][9] = 1 then error "cannot evaluate the solution further right of %1, probably a singularity", evalf[8](_val) elif _dat[4][9] = 2 then error "cannot evaluate the solution further right of %1, maxfun limit exceeded (see ?dsolve,maxfun for details)", evalf[8](_val) elif _dat[4][9] = 3 then if _dat[4][25] = 3 then error "cannot evaluate the solution past the initial point, problem may be initially singular or improperly set up" else error "cannot evaluate the solution past the initial point, problem may be complex, initially singular or improperly set up" end if elif _dat[4][9] = 4 then error "cannot evaluate the solution further right of %1, accuracy goal cannot be achieved with specified 'minstep'", evalf[8](_val) elif _dat[4][9] = 5 then error "cannot evaluate the solution further right of %1, too many step failures, tolerances may be too loose for problem", evalf[8](_val) elif 10 < _dat[4][9] then if _Env_dsolve_nowarnstop <> true then `dsolve/numeric/warning`(StringTools:-FormatMessage("cannot evaluate the solution further right of %1, event #%2 triggered a halt", evalf[8](_val), round(_dat[3][1][_dat[4][9]-10, 1]))) end if; Rounding := 'nearest'; _xout := _val else error "cannot evaluate the solution further right of %1", evalf[8](_val) end if elif _i = 2 and _xout < _val then Rounding := infinity; if _dat[4][9] = 1 then error "cannot evaluate the solution further left of %1, probably a singularity", evalf[8](_val) elif _dat[4][9] = 2 then error "cannot evaluate the solution further left of %1, maxfun limit exceeded (see ?dsolve,maxfun for details)", evalf[8](_val) elif _dat[4][9] = 3 then if _dat[4][25] = 3 then error "cannot evaluate the solution past the initial point, problem may be initially singular or improperly set up" else error "cannot evaluate the solution past the initial point, problem may be complex, initially singular or improperly set up" end if elif _dat[4][9] = 4 then error "cannot evaluate the solution further left of %1, accuracy goal cannot be achieved with specified 'minstep'", evalf[8](_val) elif _dat[4][9] = 5 then error "cannot evaluate the solution further left of %1, too many step failures, tolerances may be too loose for problem", evalf[8](_val) elif 10 < _dat[4][9] then if _Env_dsolve_nowarnstop <> true then `dsolve/numeric/warning`(StringTools:-FormatMessage("cannot evaluate the solution further left of %1, event #%2 triggered a halt", evalf[8](_val), round(_dat[3][1][_dat[4][9]-10, 1]))) end if; Rounding := 'nearest'; _xout := _val else error "cannot evaluate the solution further left of %1", evalf[8](_val) end if end if; if _EnvInFsolve = true then _dig := _dat[4][26]; _dat[4][26] := _EnvDSNumericSaveDigits; _Env_dsolve_SC_native := true; if _dat[4][25] = 1 then _i := 1; _dat[4][25] := 2 else _i := _dat[4][25] end if; _val := `dsolve/numeric/SC/IVPval`(_dat, _xout, _src); _dat[4][25] := _i; _dat[4][26] := _dig; [_xout, seq(_val[_vmap[_i]], _i = 1 .. _n-_ne)] else Digits := _dat[4][26]; _val := `dsolve/numeric/SC/IVPval`(eval(_dat, 2), _xout, _src); [_xout, seq(_val[_vmap[_i]], _i = 1 .. _n-_ne)] end if end proc, (2) = Array(0..0, {}), (3) = [t, x1(t), diff(x1(t), t)], (4) = []}); _vars := _dat[3]; _pars := map(rhs, _dat[4]); _n := nops(_vars)-1; _solnproc := _dat[1]; if not type(_xout, 'numeric') then if member(x_rkf45, ["start", 'start', "method", 'method', "left", 'left', "right", 'right', "leftdata", "rightdata", "enginedata", "eventstop", 'eventstop', "eventclear", 'eventclear', "eventstatus", 'eventstatus', "laxtol", 'laxtol', "numfun", 'numfun', NULL]) then _res := _solnproc(convert(x_rkf45, 'string')); if 1 < nops([_res]) then return _res elif type(_res, 'array') then return eval(_res, 1) elif _res <> "procname" then return _res end if elif member(x_rkf45, ["last", 'last', "initial", 'initial', "parameters", 'parameters', "initial_and_parameters", 'initial_and_parameters', NULL]) then _xout := convert(x_rkf45, 'string'); _res := _solnproc(_xout); if _xout = "parameters" then return [seq(_pars[_i] = _res[_i], _i = 1 .. nops(_pars))] elif _xout = "initial_and_parameters" then return [seq(_vars[_i+1] = [_res][1][_i+1], _i = 0 .. _n), seq(_pars[_i] = [_res][2][_i], _i = 1 .. nops(_pars))] else return [seq(_vars[_i+1] = _res[_i+1], _i = 0 .. _n)] end if elif type(_xout, `=`) and member(lhs(_xout), ["initial", 'initial', "parameters", 'parameters', "initial_and_parameters", 'initial_and_parameters', NULL]) then _xout := convert(lhs(x_rkf45), 'string') = rhs(x_rkf45); if type(rhs(_xout), 'list') then _res := _solnproc(_xout) else error "initial and/or parameter values must be specified in a list" end if; if lhs(_xout) = "initial" then return [seq(_vars[_i+1] = _res[_i+1], _i = 0 .. _n)] elif lhs(_xout) = "parameters" then return [seq(_pars[_i] = _res[_i], _i = 1 .. nops(_pars))] else return [seq(_vars[_i+1] = [_res][1][_i+1], _i = 0 .. _n), seq(_pars[_i] = [_res][2][_i], _i = 1 .. nops(_pars))] end if elif type(_xout, `=`) and member(lhs(_xout), ["eventdisable", 'eventdisable', "eventenable", 'eventenable', "eventfired", 'eventfired', "direction", 'direction', NULL]) then return _solnproc(convert(lhs(x_rkf45), 'string') = rhs(x_rkf45)) elif _xout = "solnprocedure" then return eval(_solnproc) elif _xout = "sysvars" then return _vars end if; if procname <> unknown then return ('procname')(x_rkf45) else _ndsol; _ndsol := pointto(_dat[2][0]); return ('_ndsol')(x_rkf45) end if end if; try _res := _solnproc(_xout); [seq(_vars[_i+1] = _res[_i+1], _i = 0 .. _n)] catch: error end try end proc](/view.aspx?SI=130465/430192/e4fc7fdcd6c6d9e863cc98a67040065f.gif)
|
(18) |
The output is given as a procedure. Given input t, this procedure outputs the numerical solution to the equation. For instance:
|
![[t = 2., x1(t) = HFloat(5.362531883901005), diff(x1(t), t) = HFloat(1.6948079841304207)]](/view.aspx?SI=130465/430192/cdfbc0eec870d18fee42fc2f10c5a7b0.gif)
|
(19) |
The output is a list with three items and the function x1(t) is second in this list.
To extract x1(t) from this output, we use the `op' function:
|
|
(20) |
|
|
(21) |
![MovingPoint1 := proc (t) pointplot([[x1sol(t), y1sol(t)]], symbol = solidcircle, symbolsize = 25) end proc:](/view.aspx?SI=130465/430192/3d637cc15107a4a6c9915083a03905c2.gif)
![animate(MovingPoint1, [t], t = 0 .. 5.6757, frames = 75, background = parabola);](/view.aspx?SI=130465/430192/0ad97d80babb6b266c32cb7c1569bd84.gif)
To see the point mass go back and forth, view on "continuous cycle."
For future reference, we note that the object first passes the point (3.5,0.35) at approximately t =1.46777 and again passes this point (moving in the opposite direction) at approximately t=4.22958.
|
|
(22) |
|
|
(23) |
Part 2 (in 3D) A point mass on a helix.
![f2 := [cos(t), sin(t), t];](/view.aspx?SI=130465/430192/c6860ab60d1cbbb5be3dee6f0baa2911.gif)
|
![[cos(t), sin(t), t]](/view.aspx?SI=130465/430192/54c93511c6b1f117731fef3f188bf0cb.gif)
|
(24) |
Suppose a point mass is released from rest at the point (cos(35), sin(35), 35) and travels along the helix to the xy-plane. To parameterize the motion of the point mass along the helix, we again look at the Lagrangian and solve the corresponding Euler-Lagrange equation.
Let r2(t)=[x2(t),y2(t),z2(t)] be the required parameterization.
|
|
(25) |
|
|
(26) |
|
|
(27) |
|
|
(28) |
|
|
(29) |
|
|
(30) |
Let Z denote dz/dt. The Euler-Lagrange equation d/dt (∂L/∂Z)=∂L/∂z is simply:
|
|
(31) |
We see that the height of an object falling along a helix, like an object in free fall, is a quadratic
function of time:
![z2(t) = -(1/4)*g*t^2+C[1]*t+C[2];](/view.aspx?SI=130465/430192/c4d1dc811785146d4f46a22daf70fecd.gif)
|
![z2(t) = -(1/4)*g*t^2+C[1]*t+C[2]](/view.aspx?SI=130465/430192/bc7e978d96c22b34356a6b9307b567e8.gif)
|
(32) |
The initial height is ![C[2] := 35;](/view.aspx?SI=130465/430192/96afbf6d1fd1009ce1a9d90d26c40b9e.gif)
|
|
(33) |
The initial velocity is ![C[1] := 0;](/view.aspx?SI=130465/430192/d8e2c14123f73bc54375fd005ba43d67.gif)
|
|
(34) |
And 
|
|
(35) |
|
|
(36) |
|
|
(37) |
![MovingPoint2 := proc (t) pointplot3d([[x2(t), y2(t), z2(t)]], symbol = solidcircle, symbolsize = 25, color = black) end proc:](/view.aspx?SI=130465/430192/c83bdf536911822ee9370c7554111b9d.gif)
Part 3 Motion along a piecewise defined curve.
We attach a helix to the parabola z=f1.
![h3 := [sin(t)+3.5, -cos(t)+1, t+1/4];](/view.aspx?SI=130465/430192/f7648612583c52607ff547628d6f5d7b.gif)
|
![[sin(t)+3.5, -cos(t)+1, t+1/4]](/view.aspx?SI=130465/430192/5f29805009411a46ee5b05906d3d2017.gif)
|
(38) |
|
![PLOT3D(CURVES([[3.500000000, 0., .2500000000], [3.682642477, 0.1682060356e-1, .4336734694], [3.859140641, 0.6671654885e-1, .6173469388], [4.023556880, .1480092760, .8010204082], [4.170360033, .2579639996, .9846938776], [4.294611466, .3928817100, 1.168367347], [4.392131210, .5482236121, 1.352040816], [4.459638583, .7187638173, 1.535714285], [4.494862556, .8987651471, 1.719387754], [4.496618151, 1.082172139, 1.903061223], [4.464846309, 1.262814762, 2.086734692], [4.400615873, 1.434615979, 2.270408161], [4.306087632, 1.591796189, 2.454081630], [4.184441630, 1.729067662, 2.637755099], [4.039770185, 1.841812418, 2.821428568], [3.876940220, 1.926237589, 3.005102037], [3.701429531, 1.979503009, 3.188775506], [3.519142509, 1.999816765, 3.372448975], [3.336211511, 1.986495479, 3.556122444], [3.158790554, 1.939987295, 3.739795913], [2.992848295, 1.861856803, 3.923469382], [2.843967232, 1.754732408, 4.107142851], [2.717155903, 1.622217904, 4.290816320], [2.616680394, 1.468771239, 4.474489789], [2.545920823, 1.299554543, 4.658163258], [2.507257627, 1.120260471, 4.841836727], [2.501991483, .9369206916, 5.025510196], [2.530299549, .7557029762, 5.209183665], [2.591229509, .5827037077, 5.392857134], [2.682731606, .4237427900, 5.576530603], [2.801727598, .2841678605, 5.760204072], [2.944214316, .1686743881, 5.943877541], [3.105398335, 0.8114771270e-1, 6.127551010], [3.279857230, 0.2453233728e-1, 6.311224479], [3.461721995, 0.7328714022e-3, 6.494897948], [3.644874477, 0.1054995783e-1, 6.678571417], [3.823153207, 0.5365333793e-1, 6.862244886], [3.990560673, .1285929620, 7.045918355], [4.141465086, .2328477705, 7.229591824], [4.270789839, .3629105060, 7.413265293], [4.374184292, .5144057011, 7.596938762], [4.448170130, .6822368744, 7.780612231], [4.490258381, .8607579828, 7.964285700], [4.499033144, 1.043963361, 8.147959169], [4.474199227, 1.225689758, 8.331632638], [4.416592071, 1.399823680, 8.515306107], [4.328149652, 1.560507050, 8.698979576], [4.211847280, 1.702334287, 8.882653045], [4.071597505, 1.820534150, 9.066326514], [3.912118501, 1.911130255, 9.249999983]], COLOUR(RGB, 1.00000000, 0., 0.)))](/view.aspx?SI=130465/430192/fec8a5d5a45d6298bf8e427affe7ee8a.gif)
|
(39) |
![P := spacecurve([t, 0, (t-3)^2], t = 0 .. 3.5, color = red);](/view.aspx?SI=130465/430192/94e888c81821072dbc1a792a8a1af17f.gif)
|
![PLOT3D(CURVES([[0., 0., 9.], [0.7142857143e-1, 0., 8.576530612], [.1428571429, 0., 8.163265306], [.2142857143, 0., 7.760204082], [.2857142857, 0., 7.367346939], [.3571428571, 0., 6.984693878], [.4285714285, 0., 6.612244898], [.4999999999, 0., 6.250000001], [.5714285713, 0., 5.897959184], [.6428571427, 0., 5.556122450], [.7142857141, 0., 5.224489797], [.7857142855, 0., 4.903061225], [.8571428569, 0., 4.591836736], [.9285714283, 0., 4.290816328], [.9999999997, 0., 4.000000001], [1.071428571, 0., 3.719387757], [1.142857142, 0., 3.448979595], [1.214285713, 0., 3.188775515], [1.285714284, 0., 2.938775516], [1.357142855, 0., 2.698979599], [1.428571426, 0., 2.469387763], [1.499999997, 0., 2.250000009], [1.571428568, 0., 2.040816336], [1.642857139, 0., 1.841836745], [1.714285710, 0., 1.653061236], [1.785714281, 0., 1.474489807], [1.857142852, 0., 1.306122461], [1.928571423, 0., 1.147959196], [1.999999994, 0., 1.000000012], [2.071428565, 0., .8622449099], [2.142857136, 0., .7346938893], [2.214285707, 0., .6173469502], [2.285714278, 0., .5102040927], [2.357142849, 0., .4132653166], [2.428571420, 0., .3265306220], [2.499999991, 0., .2500000090], [2.571428562, 0., .1836734775], [2.642857133, 0., .1275510274], [2.714285704, 0., 0.8163265894e-1], [2.785714275, 0., 0.4591837194e-1], [2.857142846, 0., 0.2040816645e-1], [2.928571417, 0., 0.5102042469e-2], [2.999999988, 0., 0.1440000025e-15], [3.071428559, 0., 0.5102039041e-2], [3.142857130, 0., 0.2040815959e-1], [3.214285701, 0., 0.4591836165e-1], [3.285714272, 0., 0.8163264522e-1], [3.357142843, 0., .1275510103], [3.428571414, 0., .1836734569], [3.499999985, 0., .2499999850]], COLOUR(RGB, 1.00000000, 0., 0.)))](/view.aspx?SI=130465/430192/8ff7f3433efcc3d8603bca27d94fc17d.gif)
|
(40) |
The motion of the point mass along the parabola is given by [x1sol(t),0,y1sol(t)]. When x1sol(t)=3.5, the point mass begins motion along the helix. As previously mentioned, this occurs at approximately t = 1.46777 (from line (22)).
The helix is parameterized by [sin(z3-.25)+3.5, -cos(z3-.25)+1, z3]. The initial point of motion along the helix is (3.5,0,.25).
The initial rates of change are (approximately) dx1/dt=dz3/dt = 9.26008073820564. 
As in Part 2, z3(t) is quadratic:
|
|
(41) |
|
|
(42) |
|
|
(43) |
|
|
(44) |
|
|
(45) |
So the object will first pass the point (3.5,0.25) at approximately t=1.46777, and again (in the opposite direction) at t=5.247394792. When the object returns to the parabola, it is located at the point (x1(t),y1(t)) for t = 4.22958 (from (23)).
|
|
(46) |
|
|
(47) |
|
|
(48) |
|
|
(49) |
![MovingPoint3 := proc (t) pointplot3d([[X(t), Y(t), Z(t)]], symbol = solidcircle, symbolsize = 25, color = black) end proc; -1; PH := display({H, P}); -1; animate(MovingPoint3, [t], t = 0 .. 6.6, frames = 75, background = PH)](/view.aspx?SI=130465/430192/fd3ece7bb04ff198f1874dc1520397c7.gif)
To see the point mass go back and forth, view on "continuous cycle."
Legal Notice: Maplesoft, a division of Waterloo Maple Inc. 2009. Maplesoft and Maple are trademarks of Waterloo Maple Inc. This application may contain errors and Maplesoft is not liable for any damages resulting from the use of this material. This application is intended for non-commercial, non-profit use only. Contact Maplesoft for permission if you wish to use this application in for-profit activities.
|
