08-04-01.html

08-04-01.mw

> R:=0.3:
xphug:= [ diff(theta(t),t) = ( v(t)^2 - cos(theta(t))) / v(t),
        diff(v(t),t)     = -sin(theta(t)) - R*v(t)^2 ,
        diff(x(t),t)     = v(t)*cos(theta(t)),
        diff(y(t),t)     = v(t)*sin(theta(t))];

(1)

> with(DETools):with(plots):

>
DEplot(xphug, [ theta(t), v(t), x(t), y(t) ], t=0..15,

> [seq([theta(0)=(Pi*(i/10)), v(0)=2, x(0)=0, y(0)=1 ], i=-2..2)],

> theta=-Pi/2..3*Pi, v=0..3, x=-1..6, y=-1..3,
scene=[x,y], title="path of glider",

> linecolor=[seq(COLOR(HUE,i/11),i=-2..2)]);

Plot_2d

> R:=0.3:
motorphug:= [ diff(theta(t),t) = ( v(t)^2 - cos(theta(t))) / v(t),
        diff(v(t),t)     = -sin(theta(t)) - R*v(t)^2 +k,
        diff(x(t),t)     = v(t)*cos(theta(t)),
        diff(y(t),t)     = v(t)*sin(theta(t))];

(2)

> k:=1;DEplot(motorphug, [ theta(t), v(t), x(t), y(t) ], t=0..15,

> [seq([theta(0)=0, v(0)=1+i/3, x(0)=0, y(0)=1 ], i=-5..5)],

> theta=-Pi/2..3*Pi, v=0..3, x=-1..6, y=-1..3,
scene=[theta,v], title="motor phase",

> linecolor=[seq(COLOR(HUE,i/11),i=-5..5)]);

Warning, plot may be incomplete, the following errors(s) were issued:
cannot evaluate the solution past the initial point, problem may be complex, initially singular or improperly set up

Plot_2d

> k:=0.1;DEplot(motorphug, [ theta(t), v(t), x(t), y(t) ], t=0..15,

> [seq([theta(0)=0, v(0)=1+i/3, x(0)=0, y(0)=1 ], i=-5..5)],

> theta=-Pi/2..3*Pi, v=0..3, x=-1..6, y=-1..3,
scene=[theta,v], title="motor phase, weak motor",

> linecolor=[seq(COLOR(HUE,i/11),i=-5..5)]);

Warning, plot may be incomplete, the following errors(s) were issued:
cannot evaluate the solution past the initial point, problem may be complex, initially singular or improperly set up

Plot_2d

> ?

> sol:=dsolve({op(xphug), theta(0)=0, v(0)=2, x(0)=0, y(0)=1},
numeric, range=0..1);

(3)

> sol(5);

(4)

> sol(6);

(5)

> sol(7);

(6)

> sol(8);

(7)

> sol(7.5);

(8)

> sol(7.25);

(9)

> sol(7.125);

(10)

> sol(7.1);

(11)

> sol(7.05);

(12)

> yt:= T->subs(sol(T), y(t));

(13)

> yt(7.02);

(14)

> bisector := proc( f, low, high, epsilon)
  local mid, lo, hi;
  lo := low;
  hi := high;
  mid := (lo+hi)/2;
  while ( abs(f(mid))>epsilon) do
     if ( f(mid) > 0) then
        hi:= mid;
     else
        lo:= mid;
     fi;
     mid := (lo+hi)/2;
  od;
  return(mid);
end:

let's test it for something we know the answer.

> bisector(2*x-1, 0, 1, 0.001);

Error, (in bisector) cannot determine if this expression is true or false: 0.1e-2 < abs(2*x(1/2)-1)

> bisector(2*x-1, 0, 5, 0.001);

Error, (in bisector) cannot determine if this expression is true or false: 0.1e-2 < abs(2*x(5/2)-1)

> bisector(x->2*x-1, 0, 5, 0.001);

(15)

> evalf(%);

(16)

Ramon prefers decimals

> bisector := proc( f, low, high, epsilon)
  local mid, lo, hi;
  lo := evalf(low);
  hi := evalf(high);
  mid := evalf((lo+hi)/2);
  while ( abs(f(mid))>epsilon) do
     if ( f(mid) > 0) then
        hi:= mid;
     else
        lo:= mid;
     fi;
     mid := (lo+hi)/2;
  od;
  return(mid);
end:

> bisector(x->x^2-1, 0, 5, 0.001);

(17)

> bisector(yt, 8, 7, 0.00001);

(18)

> yt(7.081542969);

(19)

> bisector(yt, 8, 7, 0.001);

(20)

Don't try this at home, kids!
it runs forever and is hard to interrupt.

> # bisector(yt, 7, 8, 0.001);

Warning, computation interrupted

> bisector := proc( f, low, high, epsilon)
  local mid, lo, hi;
  lo := evalf(low);
  hi := evalf(high);
  if (abs(f(lo))>0) then error("f not neg at lo=",lo); fi;
  if (abs(f(hi))>0) then error("f not pos at hi=",hi); fi;
  mid := evalf((lo+hi)/2);
  while ( abs(f(mid))>epsilon) do
     if ( f(mid) > 0) then
        hi:= mid;
     else
        lo:= mid;
     fi;
     mid := (lo+hi)/2;
print( [lo, hi], mid, "f(mid)=", f(mid));
  od;
  return(mid);
end:

> bisector(yt, 7, 8, 0.001);

Error, (in bisector) f not neg at lo=, 7.