1

I need to fit an implicitly defined curve. I thought I could get some data out of Solve, and then using FindFit.

Therefore, I would like to find the relation the parametric curve defined by $F(x,y)=0$:

Solve[-(1/2) + 1/2 (0.41202 BesselK[0, 0.1 Sqrt[x^2 + y^2]] + 
(0.101483 x BesselK[1, 0.1 Sqrt[x^2 + y^2]])/Sqrt[x^2 + y^2]) == 0, y]

But I can't get an output:

Solve was unable to solve the system with inexact coefficients or the system obtained by direct rationalization of inexact numbers present in the system. Since many of the methods used by Solve require exact input, providing Solve with an exact version of the system may help. >>

Edit:

In particular, I would like to fit the data coming from the curve with the expression of another curve, and not with a function $f(x)$. In particular, since this clearly looks like a cardioid, I would like it to fit to something like it.

What other strategies could I try?

usumdelphini
  • 996
  • 4
  • 14
  • For the record, "providing Solve with an exact version of the system" means using fractions instead of decimals, and the function Rationalize turns a decimal expression into one with fractions. However, Mathematica still can't solve the equation in this case. – Michael Seifert Nov 24 '15 at 18:25
  • I've plotted the contour you've defined, and it doesn't appear to be a well-defined function for $x<0$ (i.e., there are two positive values for $y$ corresponding to some negative $x$-values.) Your garden-variety fitting algorithm is going to have some trouble with that. Maybe if you wrote it in polar coordinates and solved for $\theta(r)$ instead? – Michael Seifert Nov 24 '15 at 18:32
  • Also: what's the parameter you want to adjust in your fit? Right now your curve doesn't have any free parameters for FindFit (or whatever you want to use) to adjust. – Michael Seifert Nov 24 '15 at 18:35
  • 1
    Seen this? – J. M.'s missing motivation Nov 24 '15 at 18:37

1 Answers1

4

You can get an approximation with NDSolve fairly easily for each of the four branches of the curve for which $y$ is a function of $x$. You can also get symbolic solutions in polar coordinates for each of the two branches for with $\theta$ is a function of $r$ (the reverse of the usual relationship sought).

First the polar:

eqOP = Rationalize[-(1/2) + 
     1/2 (0.41202 BesselK[0, 0.1 Sqrt[x^2 + y^2]] + 
      (0.101483 x BesselK[1, 0.1 Sqrt[x^2 + y^2]])/Sqrt[x^2 + y^2]) == 0,
   0];

tsol = Solve[
   eqOP /. {Sqrt[x^2 + y^2] -> r, 1/Sqrt[x^2 + y^2] -> 1/r} /. 
    x -> r Cos[t], t] /. C[1] -> 0

Mathematica graphics

ParametricPlot[r {Cos[t], Sin[t]} /. tsol,
  {r, 0, x /. First@NSolve[0 < x < 4 && eqOP /. y -> 0, x]}]

Mathematica graphics

wp = $MachinePrecision;
ics = {{1, y /. FindRoot[eqOP /. x -> 1, {y, 1.5}, WorkingPrecision -> wp]},
   {1, y /. FindRoot[eqOP /. x -> 1, {y, -1.5}, WorkingPrecision -> wp]},
   {x /. FindRoot[eqOP /. y -> 1/2, {x, -0.01}, WorkingPrecision -> wp], 1/2},
   {x /. FindRoot[eqOP /. y -> -1/2, {x, -0.01}, WorkingPrecision -> wp], -1/2}};
sols = Join @@ (NDSolve[{
        D[eqOP /. y -> y[x], x],
        y[First[#]] == Last[#]},
       y, {x, -1, 3},
       WorkingPrecision -> wp, PrecisionGoal -> 10, 
       AccuracyGoal -> 10, Method -> "ExplicitRungeKutta", 
       "ExtrapolationHandler" -> {Indeterminate &, 
         "WarningMessage" -> False}] & /@ ics);

(* many stiffness errors at the ends of the domains *)

Plot[y[x] /. sols // Evaluate, {x, -0.15, 2.65}]

Mathematica graphics

If you really want to fit a formula, then there is data in the solutions that can be used. For instance

Transpose@Append[
  y["Coordinates"] /. First[sols],
  y["ValuesOnGrid"] /. First[sols]]

gives the $(x,y)$ pairs for the first solution. But perhaps the InterpolatingFunction is sufficient as is.

Michael E2
  • 235,386
  • 17
  • 334
  • 747
  • Thank you very much Micheal, very very helpful! I actually wanted to fit the data coming from the curve with the expression of a curve, and not with a function $f(x)$. In particular, since this looks clearly like a cardioid, I would like it to fit to something like it. – usumdelphini Nov 25 '15 at 08:54
  • @usumdelphini You're welcome. The smooth "dimple" instead of a cusp makes it look more like a limaçon instead of a cardioid (but maybe you don't make that technical distinction), in case that helps. – Michael E2 Nov 25 '15 at 12:32