How to find the derivative value at $(\pi,0)$ for this implicit function $n$ times?

Question

I am trying to take the implicit derivative at $\sin(x+y)+\sin(x)=y$ and substitute $x=\pi$ and $y=0$ at least 6-7 times since I need to find the Taylor series for this function.

Since I barely anything about coding I referred to the documentation in Mathematica and other web pages which showed to do this...

D[Sin[x]+Sin[x+y[x]] == y[x],{ x,1}]

Then I took the output of that and solved for $y[x]$

Solve[Cos[x] + Cos[x + y[x]] (1 + Derivative[1][y][x]) == 
 Derivative[1][y][x], Derivative[1][y][x]]

Then took that output and finally substituted $x=\pi$ and $y[x]=0$

    {Derivative[1][y][x] -> (-Cos[x] - Cos[x + y[x]])/(-1 + 
        Cos[x + y[x]])}/.{y[x]->0,x->Pi}
(*-1*)

Unfortunately this process was too tedious so I took some programming from here

D[Sin[x] + Sin[x+y[x]] == y, {x, 1}]
sol1 = Solve[%, y'[x]]
D[Sin[x] + Sin[x+y[x]] == y, {x, 2}]
sol2 = Solve[%, y''[x]]
sol2 /. sol1 // Simplify

And ended up with the second derivative but was too long.

Is there a more elegant way of doing this and how? (One thing you can do is solve for the first implicit derivative, substitute $x$ and $y$, then substitute the first into the second derivative and so on with the next derivatives).

Answer 1

You can substitute a Taylor expansion in your equation

y=Sum[a[n](x-π)^n, {n,1,6}]+ O[x,π]^7

SolveAlways[ Sin[x+y]+Sin[x]-y ==0, x]

(* {{a[5]->-(1/240),a[6]->0,a[4]->0,a[3]->1/12,a[2]->0,a[1]->-1}} *)

Answer 2

Last[Nest[Function[d, {d, Append[Last[#], Solve[0 == d /. x -> π /. Last[#]][[1, 1]]]}][
                   D[First[#], x]] &, {Sin[x] + Sin[x + y[x]] - y[x], {y[π] -> 0}}, 15]]

$$\left\{y(\pi )\to 0,y'(\pi )\to -1,y''(\pi )\to 0,y^{(3)}(\pi )\to \frac{1}{2},y^{(4)}(\pi )\to 0,y^{(5)}(\pi )\to -\frac{1}{2},y^{(6)}(\pi )\to 0,y^{(7)}(\pi )\to \frac{1}{2},y^{(8)}(\pi )\to 0,y^{(9)}(\pi )\to 17,y^{(10)}(\pi )\to 0,y^{(11)}(\pi )\to -\frac{1153}{4},y^{(12)}(\pi )\to 0,y^{(13)}(\pi )\to \frac{13297}{4},y^{(14)}(\pi )\to 0,y^{(15)}(\pi )\to \frac{108109}{2}\right\}$$

Answer 3

You need to find the Taylorpolynom for your function?

The code for the Taylor formula:

taylor = (vars - point).# &;
init := D[f[vars], {vars, j}] /. Thread[vars -> point];
taylorPolynom[m_] := Sum[1/j! Nest[taylor, init, j], {j, 0, m}]

vars = {x, y};
point = {Pi, 0};
f[vars_] = Sin[x + y] + Sin[x] - y;

taylorPolynom[2]
-2 (-Pi + x) - 2 y

or

taylorPolynom[6] // FullSimplify
1/120 (2 (120 + (-20 + (Pi - x)^2) (Pi - x)^2) (Pi - x) - 
   5 (48 + (-12 + (Pi - x)^2) (Pi - x)^2) y + 
   10 (-6 + (Pi - x)^2) (Pi - x) y^2 - 
   10 (-2 + (Pi - x)^2) y^3 + 5 (Pi - x) y^4 - y^5)

Answer 4

You can use the new in M12 function AsymptoticSolve to do this:

AsymptoticSolve[Sin[x + y] + Sin[x] == y, {y, 0}, {x, π, 10}]

{{y -> π - x + 1/12 (-π + x)^3 - 1/240 (-π + x)^5 + (-π + x)^7/10080 + (17 (-π + x)^9)/362880}}

Answer 5

The SolveAlways form is very nice, while some of the other solns didn't immediately work V7. But what is SolveAlways doing? Probably something like the following -- Compute the derivatives by brute force, and then match at the expansion center. The derivative computation becomes increasingly time consuming as the \ number of terms increases, but 9 terms is pretty fast.
The Taylor series form is used here, so the y derivatives are obtained, as shown by Coolwater. Removing the n! will match the other \ responses. (Note, subscripts are inert labels.)

nterm = 9;
Clear[yb, dyb, aList]
yb = Sum[Subscript[a, n] (x - Pi)^n/n!, {n, 1, nterm + 1}] + O[x, Pi]^(nterm + 2);
dyb = Table[   D[Sin[x + yb] + Sin[x] - yb, {x, jj}] /. x -> Pi, {jj, 1, nterm}];
aList = Table[Subscript[a, n], {n, 1, nterm}];
(* match the representation at Pi by adjusting the coefficients *)
vDerivs = Solve[dyb == 0, aList][[1]];
Print[vDerivs];
(* {Subscript[a, 8]->0,Subscript[a, 9]->17,Subscript[a, 7]->1/2,Subscript[a, 6]->0,Subscript[a, 5]->-(1/2),Subscript[a, 4]->0,Subscript[a, 3]->1/2,Subscript[a, 2]->0,Subscript[a, 1]->-1} *)

Answer 6

I'm surprised that nobody followed up on QuantumDot's suggestion to use Dt[]:

With[{n = 9},
     (Table[Dt[y, {x, k}], {k, n}] /. 
      First @ Solve[Table[Dt[Sin[x] + Sin[x + y] - y, {x, k}] == 0, {k, n}], 
                    Table[Dt[y, {x, k}], {k, n}]]) /. {x -> π, y -> 0}]
   {-1, 0, 1/2, 0, -1/2, 0, 1/2, 0, 17}