0

The code is supposed to move the center of a circle of radius "d" incrementally on a function (the one in the first line), then solve for the two intercepts and repeat. I noticed that it gave two complex coordinates for every point as well (not really sure why, Id be glad if someone could explain), but I ignored it. However, as soon as the program reaches around halfway through the domain (x0 = 320), it outputs only complex answers, but logically, the circle should still intercept the function twice. What's going on?

y[x_] := (x Tan[Theta]) - ((10*(x^2))/(2*(v^2)*((Cos[Theta])^2)))
a = Table[sol := Solve[y1 == y[x1] &&
                       d == EuclideanDistance[{x0, y0}, {x1, y1}], {x1, y1}];
z = N[sol //. {x0 -> i, v -> 80, Theta -> Pi/4, d -> 0.8255, y0 -> y[x0]}], {i, 0, 640, 2}]
J. M.'s missing motivation
  • 124,525
  • 11
  • 401
  • 574
Agent Wowza
  • 3
  • 2
  • 1
    Can you edit your question to add some context? In particular, why do you expect real solutions to exist for all $x_0, y_0$? Right now it just looks like the system$$y_1 = x_1 - \alpha x_1^2 \(x_1-x_0)^2 + (y_1 - y_0)^2 = d^2,$$which can be reduced to a fourth-order real polynomial in $x_1$. It's not surprising that depending on the values of the coefficients of this polynomial, there can be four, two, or zero real roots. – Michael Seifert Aug 29 '17 at 14:58
  • You can call Solve[eq,vars,domain] with Reals as a third argument to limit the solution domain to strictly real solutions. Also you can make solving much faster by applying your replacement rules before solving. You also only need to solve once and then you can reuse the general solution obtained by Solve. – Thies Heidecke Aug 29 '17 at 14:58
  • @MichaelSeifert As he explained, he's basically moving the center of a circle along a curve. Geometrically, you'd expect the circle to intersect the curve at least twice, provided the curve does not have discontinuities – Lukas Lang Aug 29 '17 at 15:07
  • @MichaelSeifert I edited the question to show that I'm only looking at the domain 0<x0<640, but yes as Mathe172 says, geometrically it makes sense. – Agent Wowza Aug 29 '17 at 15:09
  • @ThiesHeidecke Edit- Sorry spoke too soon, I tried the Reals argument but evaluation is taking way too long now, even though after abortion, the program seems to have moved on. – Agent Wowza Aug 29 '17 at 15:10
  • @AgentWowza Look at my answer to see an example how to put the Reals argument to use without the long waiting time. The intermediate answers you get, when you paste the code into Mathematica and evaluate it are also helpful for understanding. – Thies Heidecke Aug 29 '17 at 15:25

2 Answers2

2

Solution

I tried to take your code and simplify it a bit to make things easier to understand and faster:

y[x_] = (x Tan[theta]) - ((10*(x^2))/(2*(v^2)*((Cos[theta])^2)))
initialconditions = {x0 -> i, v -> 80, theta -> Pi/4, d -> Rationalize[0.8255], y0 -> y[x0]}
eqs = y1 == y[x1] && d^2 == SquaredEuclideanDistance[{x0, y0}, {x1, y1}]
sol = Solve[eqs //. initialconditions, {x1, y1}, Reals]
a = Transpose@Table[ N[{x1, y1} /. sol], {i, 320, 640, 1}];
ListPlot[a]

Plot of two near coincident trajectories

Explanation

I'll go through every line to explain what and why i changed it:

y[x_] = (x Tan[theta]) - ((10*(x^2))/(2*(v^2)*((Cos[theta])^2)))

For the definition of y a = (Set) instead of := (SetDelayed) suffices, since the right hand expression will always evaluate to the same thing, we just want classic function behaviour here. Also i put theta in lowercase since upper case symbols are usually reserved for functions.

initialconditions = {x0 -> i, v -> 80, theta -> Pi/4, d -> Rationalize[0.8255], y0 -> y[x0]}

We can put the replacement rules in a named expression for better readability and put Rationalize around the value of d in case we are interested in exact solutions.

eqs = y1 == y[x1] && d^2 == SquaredEuclideanDistance[{x0, y0}, {x1, y1}]

Distance constraints are usually better (read easier to solve for mathematica and give simple solutions) formulated in terms of squared distance, since we can get rid of a square root this way without changing the meaning of our constraint.

sol = Solve[eqs //. initialconditions, {x1, y1}, Reals]

We can now compute the general solution outside the loop (which means we only need to compute it once), but with the initial conditions already substituted (which means Solve will give us a result much quicker). Also we used the Reals domain argument to constrain the solution set, which gives us two real solution for each i.

a = Transpose@Table[ N[{x1, y1} /. sol], {i, 320, 640, 1}];

We let Table take our general solution and let it go through every i value, which now all reduce to numerical values. The Transpose is there to bring it in a form of two lists of points, suitable to be displayed by ListPlot.

Thies Heidecke
  • 8,814
  • 34
  • 44
1

It seems the general solution returned by Solve doesn't return all/any solutions when $y'(x)=0$. You can circumvent this by inserting the values first and solving afterwards, this ensures that Solve sees the special case and finds all solutions:

y[x_] := (x Tan[Theta]) - ((10*(x^2))/(2*(v^2)*((Cos[Theta])^2))) 
a = Table[
  sol = Solve[
    (y1 == y[x1] && d == EuclideanDistance[{x0, y0}, {x1, y1}])
     //. {x0 -> i, v -> 80, Theta -> Pi/4, d -> 0.8255, y0 -> y[x0]},
    {x1, y1}
    ];
  z = N[sol],
  {i, 320, 640, 1}
  ]
Lukas Lang
  • 33,963
  • 1
  • 51
  • 97