1
$\begingroup$

I have the following piecewise continuous function:

myfun[a_, l1_, l2_] = 
 Piecewise[{{a^2, a < l1}, 
            {a - l1 + l1^2, a >= l1 && a < l2}, 
            {l2 - l1 + l1^2 + Log[a - l2 + 1], a >= l2}}]

The function is defined by a value $a$ and two parameters $l_1$ and $l_2$ and if I plot it, appears as:

Plot[Evaluate[myfun[a, l1, l2] /. {l1 -> 1.5, l2 -> 3.2}], {a, 0, 5}, PlotRange -> All]

enter image description here

Now, I would like to invert the numerical solution given by $myfun==2.0$ with respect to a, which is a function depending on $l_1$ and $l_2$. I tried to do it "in a delayed way" with

myInvFun[l1_, l2_] := x /. FindRoot[myfun[x, l1, l2] == 2.0, {l, 0}]

But this does not help me, since if I try to plot the following

Plot[Evaluate[myInvFun[l1, l2] /. {l2 -> 3.5}], {l1, 0, 3.5}]

errors appear:

enter image description here

How can I make the inverse function usable as an explicit function even if I use only a definition from numerical solution?

P.S. I had a look to answers like this one and this one but either way I could not find a way to solve my problem

$\endgroup$

closed as off-topic by Anton Antonov, user9660, MarcoB, ubpdqn, Öskå May 7 '16 at 6:46

This question appears to be off-topic. The users who voted to close gave this specific reason:

  • "This question arises due to a simple mistake such as a trivial syntax error, incorrect capitalization, spelling mistake, or other typographical error and is unlikely to help any future visitors, or else it is easily found in the documentation." – Anton Antonov, Community, MarcoB, ubpdqn, Öskå
If this question can be reworded to fit the rules in the help center, please edit the question.

5
$\begingroup$

This is simple enough to do analytically.. assuming l2>l1 you can readily find the transition points and invert each piece of the Piecewise expression:

myfuninv[x_, l1_, l2_ /; l2 > l1] = Piecewise[{
   {Sqrt[x], x <= l1^2},
   {x + l1 - l1^2, l1^2 < x <= l2 + l1^2 - l1},
   {Exp[x + l1 - l2 - l1^2] - 1 + l2,  l2 + l1^2 - l1 < x}}]
Plot[ myfuninv[2, l1, 3.5] , {l1, 0, 3.5}]

enter image description here

$\endgroup$
  • 2
    $\begingroup$ ...in the words of Field Marshal the Viscount Sir William Slim: “Thinking is the cheapest and one of the most effective long-range weapons.” (+1) $\endgroup$ – Michael E2 May 6 '16 at 16:21
4
$\begingroup$

You can see what you get with the body of your function myInvFun over parameter instances:

In[215]:= FindRoot[myfun[x, 2, 5] == 2.0, {x, 0}]

Out[215]= {x -> 0.}

FindRoot::jsing: Encountered a singular Jacobian at the point {x} = {0.}. 
Try perturbing the initial point(s). >>

So, if you change your definition as advised by the message you'll get results you might expect:

myInvFun[l1_, l2_] := x /. FindRoot[myfun[x, l1, l2] == 2.0, {x, 0.1}]

Plot[Evaluate[myInvFun[l1, 3.5]], {l1, 0, 3.5}]

enter image description here

$\endgroup$
2
$\begingroup$

General Inversion

You can use FindRoot to do a general inversion of a Piecewise function.

The strategy will be to extract the smooth continuous function from the piecewise function and use that as input to FindRoot.

Below is a copy of your function:

myfun[a_, l1_, l2_] = Piecewise[{
   {a^2, a < l1},
   {a - l1 + l1^2, 
    a >= l1 && a < l2}, {l2 - l1 + l1^2 + Log[a - l2 + 1], a >= l2}
   }]

which has three regions that contain three different functions.

Part can be used to access the function that is valid in each region

myfun[a, l1, l2][[1, 1, 1]]
(* a^2 *)

myfun[a, l1, l2][[1, 2, 1]]
(* a - l1 + l1^2 *)

myfun[a, l1, l2][[1, 3, 1]]
(* -l1 + l1^2 + l2 + Log[1 + a - l2] *)

Now we will define the inverse function which takes the limits as well as myfun as input arguments. myfun will be used to set the limits of the inverse Piecewise function as well as being an input to FindRoot.

myInvFun[y_, l1_, l2_, myfun_] := Piecewise[{
   {
    Module[{a, fun},
     fun[a_] = myfun[a, l1, l2][[1, 1, 1]];
     FindRoot[fun[a] == y, {a, l1/2}][[1, 2]]
     ],
    y < myfun[l1, l1, l2]
    },
   {
    Module[{a, fun},
     fun[a_] = myfun[a, l1, l2][[1, 2, 1]];
     FindRoot[fun[a] == y, {a, (l1 + l2)/2}][[1, 2]]
     ],
    y >= myfun[l1, l1, l2] && y < myfun[l2, l1, l2]
    },
   {
    Module[{a, fun},
     fun[a_] = myfun[a, l1, l2][[1, 3, 1]];
     FindRoot[
       l2 - l1 + l1^2 + Log[a - l2 + 1] == y, {a, l2 + 0.5}][[1, 2]]
     ],
    y >= myfun[l2, l1, l2]
    }
   }]

Now myInvFun can be treated like any ordinary function. For example:

myInvFun[2, 1.5, 3.2, myfun]
(* 1.41421 *)

We can validate that the inverse matches myfun by plotting myfun and using ParametricPlot for myInvFun.

Show[
 Plot[myfun[a, 1.5, 3.2], {a, 0, 5}, PlotRange -> All, 
  PlotStyle -> Black],
 ParametricPlot[{myInvFun[y, 1.5, 3.2, myfun], y}, {y, 0, 5}, 
  PlotStyle -> {Red, Dashed}]
 ]

Mathematica graphics

They do indeed match.

Specific Inversion

Clearly you can simply replace the argument a in myFunInv with 2 if you want to plot the inverse as function of l1.

Plot[myInvFun[2, l1, 3.5, myfun], {l1, 0, 3.5}]

Mathematica graphics

$\endgroup$
2
$\begingroup$

The following uses NDSolve to construct interpolations along lines in the domain. We then construct an interpolation between the solutions to NDSolve, which represents a as a function of l1 and l2.

{dadl1, dadl2} = 
  grad = PiecewiseExpand /@ (-D[myfun[a, l1, l2], {{l1, l2}}]/
        D[myfun[a, l1, l2], a]) /. ComplexInfinity -> 0 // Quiet;

sols = Table[
   Block[{l1a = 0, l1b = 3.5, l2},
    Module[{flag = True, res, dl = 0, n = 0},
     While[flag,
      Check[
       l2 = l2b + (-1)^n dl;
       res = 
        NDSolveValue[{a'[l1] == (dadl1 /. a -> a[l1]), a[0] == 2}, 
         a, {l1, l1a, l1b}, PrecisionGoal -> 10]; flag = False,
       flag = True;
       n += 1;
       dl += 0.001]
      ];
     {l2, res}]
    ],
   {l2b, 0., 3.5, 0.1}];

ClearAll[if];
if[sols : {{_?NumericQ, _InterpolatingFunction} ..}] := 
  if[sols[[All, 1]], sols[[All, 2]], 
   Nearest[sols[[All, 1]] -> Automatic]];
if[l2data_, ifns : {__InterpolatingFunction}, nf_NearestFunction][
  l1in_?NumericQ, l2in_?NumericQ] := Module[{i2, $l2, $a},
  {i2[1], i2[2]} = nf[l2in] /. {{1} -> {1, 2},
      {Length@l2data} -> {Length@l2data - 1, Length@l2data},
      {n_ /; l2in <= l2data[[n]]} :> {n - 1, n}, {n_} :> {n, n + 1}};
  {$l2[1], $l2[2]} = l2data[[{i2[1], i2[2]}]];
  $a[1] = ifns[[i2[1]]][l1in];
  $a[2] = ifns[[i2[2]]][l1in];
  InterpolatingPolynomial[{{{$l2[1]}, $a[1]}, {{$l2[2]}, $a[2]}}, l2in]
  ]

afn = if[sols];

Plot3D[afn[l1, l2], {l1, 0, 3.5}, {l2, 0, 3.5}]

Mathematica graphics

$\endgroup$

Not the answer you're looking for? Browse other questions tagged or ask your own question.