Message Boards Message Boards

0
|
5423 Views
|
4 Replies
|
0 Total Likes
View groups...
Share
Share this post:

Efficient way to solve parametric equations?

Posted 10 years ago
I'm trying to calculate the intersection point of a curve given by a parametric equation and a vector. (i.e. a circle and a ray vector bouncing within the circle).
So far I couldn't find an efficient way to calculate the intersection points of the vector and the curve. I tried "FindRoot", "Solve", "Reduce" etc. but all have their own drawbacks. (e.g. The problem with FindRoot was that I could not find an easy way to generate a good starting point for the parameter t. But without it, the obvious solution with a vector of length 0 at r[t0] was returned...)
A general problem of my other solutions is, that they are all extremely slow.
For further processing I need the intersection point given by the parameter t and not only x,y coordinates. I found out that Solve is MUCH faster when the curve is given in implicit form rather than in parametric form. Unfortunatly that doesn't help me (at least I don't know how). Is there any way to speed up the calculation?

r[t_] := {Cos[t], Sin[t]}(*parametric curve*)

uT[t_] :=
Simplify[r'[t]/Norm[r'[t]], t \[Element] Reals](*tangent vector at t*)

uN[t_] :=
Simplify[uT'[t]/Norm[uT[t]], t \[Element] Reals](*normal vector at t*)
 t0 = 0;(*start parameter t*)
 a0 = 45.1/180*\[Pi]; (*start parameter a; a=angle between tangent vector \
 and ray*)
 T = {t0};
 A = {a0};
 trace[t0_, a0_, n_] := Module[{t = t0, a = a0},
    Do[
     rr = RotationMatrix[a].uT[t];
     (*res2=Solve[{xx^2+yy^2==1,yy==rr[[2]]/rr[[1]]*xx+(r0[[2]]*
    rr[[1]]-rr[[2]]*r0[[1]])/rr[[1]]},{xx,yy},Reals];*)
   
    res = FindInstance[
      Reduce[{r[tt] == r[t] + k*rr, Abs[k] > 10^(-10), tt < 2 \[Pi],
        tt >= 0}, {tt, k}], {k, tt}];
    AppendTo[T, tt /. res[[1]]];
    t = tt /. res[[1]];
    rr = 2*(uT[t].rr)*uT[t] - rr;
    a = ArcCos[rr.uT[t]];
    AppendTo[A, a];
    , {i, n}]
   ];
trace[t0, a0, 10]
Reduce::ratnz: Reduce was unable to solve the system with inexact coefficients. The answer was obtained by solving a corresponding exact system and numericizing the result. >>

Reduce::ratnz: Reduce was unable to solve the system with inexact coefficients. The answer was obtained by solving a corresponding exact system and numericizing the result. >>

Reduce::ratnz: Reduce was unable to solve the system with inexact coefficients. The answer was obtained by solving a corresponding exact system and numericizing the result. >>

General::stop: Further output of Reduce::ratnz will be suppressed during this calculation. >>
Show[
ParametricPlot[{r[t]}, {t, 0, 2 \[Pi]}],
ListLinePlot[r /@ T]]
POSTED BY: Julian Stuermer
4 Replies
Posted 10 years ago
Here is code that should be adaptable for the purpose at hand. It returns a list {total time, x intersection, y intersection, new x direction, new y direction}. Hence it is suitable for nesting.
 solveBilliard[{t0_, x0_, y0_, dx_, dy_}] :=
  Catch[Module[{x, y, t, polys, soln, unitv, norm, perp},
    polys = Flatten[{x^2 + y^2 - 1,
       Thread[{x, y} - ({x0, y0} - t*{dx, dy})]}];
    soln = {t, x, y} /. NSolve[polys == 0, {t, x, y}];
    soln = Select[soln, First[#] != 0 &];
    If[Length != 1, Throw[$Failed], soln = First;
     soln[[1]] += t0];
    unitv = -{dx, dy};
   unitv = unitv/Norm;
   norm = -Rest/Norm[Rest];
   perp = unitv - Projection[unitv, norm];
   unitv = norm - perp;
   unitv = unitv/Norm;
   Join[soln, unitv]
   ]]

Could use as follows.
hitsM = NestList[solveBilliard, {0., 1., 0., -Sqrt[2]/2., Sqrt[2]/2.}, 100];
ptsM = hitsM[[All, {2, 3}]];

Show[ListPlot[ptsM, Joined -> True], Graphics[Circle[{0, 0}, 1]],  PlotRange -> All, AspectRatio -> Automatic]

One caveat is that it might require modification to use fairly high precision to correctly follow a large number of bounces. See the SIAM 100 Digit Challenge problem 2 for example.
POSTED BY: Updating Name
Posted 10 years ago
@sean:  Rationlize the fraction 45.1/180 to 451/1800 lead to an more or less infinite calculation time... I aborted the calculation after 10min.
@Daniel: Yes the curve is a simple unit circle here. This was just a test, as here calculations and results are easy to interpret.
You are right: The vector is originating on the circle and it is the rotated tangent vector. At the next intersection point the vector gets reflected etc.
It's supposed to be a billiard system ( http://mathworld.wolfram.com/Billiards.html ).

I guess the problem is related to the periodicity of the function Cos/Sin.
Interestingly the solution is found much faster when I do not restrict parameter tt to be between 0 and 2 pi. But then the solution is t=208.919 etc.
Using Mod[t, 2 pi] afterwards works. It's still not very efficient, but at least it is working. Any ideas how I can speed up it further ?
POSTED BY: Julian Stuermer
Your parametric curve is the unit circle. I was unable to discern what was the segment with which you wish to intersect that curve. Is it the rotated tangent vector originating on the circle?
POSTED BY: Daniel Lichtblau
The first thing to address is the warning message:
Reduce::ratnz: Reduce was unable to solve the system with inexact coefficients. The answer was obtained by solving a corresponding exact system and numericizing the result. >>
This message means you have given Reduce (or something related) an equation with floating point numbers. This can often be problematic (as floating point numbers tend to be).

Try defining a0 as:
a0 = 451/1800*?

What error messages do you get from Resolve then?
POSTED BY: Sean Clarke
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract