# Implicit function theorem to higher order

Consider

f[x_, y_] := Csc[0.482 y] Sin[0.963 x - 0.482 y] +
3.247 Csc[0.333 y] Sin[0.667 x - 0.333 y] +
5.049 Csc[0.119 y] Sin[0.238 x - 0.119 y]


I'd like an expression of an approximation of the root curve of $f$ (in the neighbourhood of a chosen known root). This is closely related to the implicit function theorem. The theorem gives an easy way to plot a linear approximation around a known root (point), corresponding to the red dot in the $(x,y)$ graph:

Question How can I get an better approximation (to a higher order) of the roots (blue curves)? Of course I am not looking for an interpolation of the blue curves, my goal is precisely to avoid solving many equations (each of the blue dot is the result of a NSolve). Instead, I'd like a parametrisation $y=\phi(x)$ where $\phi(x)$ could be a Taylor expansion at order $n$, and such that $f(x,\phi(x))=o(n)$.

Full code for figure:

f[x_, y_] := Csc[0.482 y] Sin[0.963 x - 0.482 y] +
3.247 Csc[0.333 y] Sin[0.667 x - 0.333 y] +
5.049 Csc[0.119 y] Sin[0.238 x - 0.119 y]
sols = Table[
Map[{#, y} &,
NSolve[f[x, y] == 0 && 0 < x <= y/2, x][[All, 1, 2]]], {y, 0,
20, 0.05}] // Flatten // Partition[#, 2] &; // Quiet
point = sols[[255]];
plot = Show[ListPlot[sols],
Graphics[{Red, PointSize[0.03], Point[point]}]]

{x0, y0} = point;
phi[x_] = y0 - (x - x0)*(D[f[x, y], x] /. x -> x0 /.
y -> y0)/(D[f[x, y], y] /. x -> x0 /. y -> y0) // Simplify
Show[plot, Plot[phi[y], {y, -10, 25}, PlotStyle -> Red]]

• It is not clear what you want for a result. A parametric InterpolatingFunction object? A function of the form y=... with the ... a series in x? – Daniel Lichtblau Sep 15 '15 at 0:59
• @DanielLichtblau The second one ($y(x)$). I edited. – anderstood Sep 15 '15 at 1:43

A first-year calculus approach to finding Taylor series:

order = 10; (* derivative order *)
step[x0_][{eqn_, coeffs_}] :=
{#, {coeffs, Solve[# /. x -> x0 /. Flatten@coeffs]}} &[D[eqn, x]];
derivatives =
Flatten@ Last@ Nest[step[x0], {f[x, y[x]] == 0, y[x0] -> y0}, order];
y1 = Normal@Series[y[x], {x, x0, Length@derivatives - 1}] /. derivatives
(*
11.4 + 1.12299 (-2.26109 + x) + 0.190117 (-2.26109 + x)^2 +
...
0.0108529 (-2.26109 + x)^9 + 0.00907805 (-2.26109 + x)^10
*)

Show[
plot,
Plot[y1, {x, 0, 3.3}, PlotStyle -> {Red, Opacity[0.5]}]
]


Notes: The function Nest applies step iteratively to the equation f[x, y[x]] == 0 to calculate the value of the next derivative of the point {x0, y0} in question. The function step differentiates the equation one more time at each step and solves for the next derivative. The derivative will be the only unknown in the equation, so we do not have to specify which unknown to solve for. Solve will figure it out for us. The function step returns the differentiated equation and the solutions for all the derivatives up to the current order; the derivative values have the form of a nested (linked) list of rules.

The divergence of the red curve on the left occurs outside the radius of convergence of the Taylor series.

• Oh, well, just noticed Daniel beat me to it. – Michael E2 Sep 15 '15 at 16:40
• Yeah, but with messier code. – Daniel Lichtblau Sep 15 '15 at 16:44

Write y as an explicit function of x. Then one can solve for successive derivatives to set up a Taylor approximation. Below is some slightly messy code for this.

taylor[func_, x_, y_, pt_, n_] := Module[
{f = func[x, y] /. y -> y[x], deriv, var = y[x], sol, newsol},
deriv = f;
sol = {y[x] -> pt[[2]]};
pt[[2]] + Sum[
deriv = D[deriv, x];
var = D[var, x];
newsol =
Solve[(deriv /. sol /. x -> pt[[1]]) == 0,
var /. x -> pt[[1]]] /. sol;
sol = Join[sol, newsol[[1]]];
1/j!*(x - pt[[1]])^j*(var /. x -> pt[[1]]) /. sol
, {j, n}]
]


Illustrated on the example in question:

plot = Show[ListPlot[sols],
Graphics[{Red, PointSize[0.03], Point[point]}],
Plot[Evaluate[taylor[f, x, y, point, 4]], {x, point[[1]] - 2,
point[[1]] + 2}, ColorFunction -> (Green &)]]


• I accepted Michael's answer because of the explanations. But your answer is (also) what I was looking for, thank you very much. – anderstood Sep 15 '15 at 16:58

Here I'll use a Chebyshev series instead of a Taylor series (see About multi-root search in Mathematica for transcendental equations).

First we approximate the curve of interest. The interpolation will be used to seed FindRoot below to get more precise values of y for a given x.

yIF = NDSolveValue[{f[x, y[x]] == 0, u'[x] == 1, u[x0] == x0,
y[x0] == y0}, y, {x, 0, 10}];


NDSolveValue::ndsz: At x == 3.275164228822655, step size is effectively zero; singularity or stiff system suspected. >>

The Chebyshev proxy approximates an analytic function very well; that means that to get take advantage of its strengths, we need to stay a little away from the singularity (at the vertical tangent where the curve loops back).

Module[{x1, x2},
{x1, x2} = First[yIF["Domain"]];
domain = {x1, x1 + 0.999 (x2 - x1)}];

r = 10;
yFR = y /. FindRoot[f[#, y] == 0, {y, yIF[#]}] &;
n = 32;  (* use n = 64 for greater accuracy *)
cnodes = Rescale[N[Cos[Pi Range[0, n]/n], 30], {-1, 1}, domain];
cc = Sqrt[2/n] FourierDCT[yFR /@ cnodes, 1];
cc[[{1, -1}]] /= 2;


To get a numerically stable result, we need to evaluate the Chebyshev polynomials individually before summing with the coefficients. (The expanded polynomial has alternating coefficients.) Note that the result yT is a polynomial.

yT[x_?NumericQ] :=
cc.Table[ChebyshevT[n - 1, Rescale[x, domain, {-1, 1}]], {n, Length@cc}];

Show[
plot,
Plot[yT[x], Evaluate@Flatten[{x, domain}],
PlotStyle -> {Red, Opacity[0.5]}]
]


• Why do you approximate the function of interest? Is it just to give good initial conditions to FindRoot? – anderstood Sep 15 '15 at 17:12
• @anderstood Yes it is just to give good initial starting points to FindRoot. Since there are many possible roots, it's important to start close to the desired value. -- Just noticed there was a typo that obscured the reason. – Michael E2 Sep 15 '15 at 17:19
• The accuracy of the fit (considering the short calculation time) and the length of the domain of "validity" are impressive. However, I will stick with the Taylor series because it still stands if I replace the values by parameters name, so it gives a closed form approximation in a more general frame. Thank you very much for taking the time to share this informative technique. – anderstood Sep 15 '15 at 18:30

The function is easier to work with, if singularities are eliminated.

f[x_, y_] == (Csc[0.482 y] Sin[0.963 x - 0.482 y] +
3.247 Csc[0.333 y] Sin[0.667 x - 0.333 y] +
5.049 Csc[0.119 y] Sin[0.238 x - 0.119 y])
Sin[0.482 y] Sin[0.333 y] Sin[0.119 y] // Simplify


Next, ContourPlot quickly finds all the zero-curves.

plt = ContourPlot[f == 0, {x, 0, 10}, {y, 0, 20}, PlotPoints -> 100]


Points on a segment of any of the curves shown in this plot or that in the Question can be extracted in a variety of ways. Then, Fit or a related Mathematica function can be used to fit these points to an analytical function y[x].

A simpler approximation can be obtained based on my answer to a related question. From it, y can be approximated locally by

y0 - (f[x, y0]/Derivative[0, 1][f][x, y0]);


With y0 = 11.4, as specified in the question, this becomes

% /. y0 -> 11.4 // FullSimplify
(* 33.4981 + (-64.2294 Cos[5.4948 - 0.963 x] -
167.835 Cos[3.7962 - 0.667 x] + 58.1109 Cos[1.3566 - 0.238 x] -
67.2963 Sin[5.4948 - 0.963 x] +
205.631 Sin[3.7962 - 0.667 x])/(2.90656 Cos[5.4948 - 0.963 x] +
7.595 Cos[3.7962 - 0.667 x] - 2.62968 Cos[1.3566 - 0.238 x] +
2.77246 Sin[5.4948 - 0.963 x] - 10.3375 Sin[3.7962 - 0.667 x] +
1. Sin[1.3566 - 0.238 x]) *)


which is a reasonable local approximation to the curve.

Show[plt, Plot[%, {x, 1.5, 3.5}, PlotStyle -> Red]]


• My question might be not clear enough, but the problem with such a solution is that I need to calculate the approximation if I change the values of the numbers in $f$. I am looking for something more in the spirit of Taylor expansion which would give a formula independently of these values. – anderstood Sep 15 '15 at 4:32
• @anderstood Is the addendum to my answer closer to what you had in mind? Converting the expression for the curve to a Taylor series is straightforward but slightly less accurate. Also, choosing y0 = 11 would give a wider range of validity for the approximation. Incidentally, expanding both the numerator and the denominator of the expression would give a result similar to a Pade approximation. – bbgodfrey Sep 15 '15 at 15:18
• That is exactly in the spirit of what I am looking for, indeed. Using Series[f[x, y], {y, y0, 3}] // Normal; ysol = y /. Solve[% == 0, {y}], I managed to extend it to the third order by expanding $f(x,y)$ with Taylor to the third order in $y$, and then solve the third order polynomial. But it is still limited as findings roots of polynomials, for which there is in general no analytical solutions if $n>5$. – anderstood Sep 15 '15 at 16:22
• @anderstood At some point, analytical solutions simply do not exist, and I think you have reached that point. An alternative is to determine the solution numerically and fit it to an analytical solution. – bbgodfrey Sep 16 '15 at 0:01
• I indeed reached that point, but for some reasons it is important to me to get something analytical, even if it is an approximation. This allows me to give an analytical expression of a surface which would be very costly to build numerically. I'm trying to do the Padé approximation you mentioned (hadn't heard about it before), but I don't think it is as easy as Taylor (see Michael's answer for example) because it does not make sense to identify successive derivatives of $x\mapsto y(x)$ in $x_0$. – anderstood Sep 16 '15 at 0:33