# About two point Taylor series expansion

From a comment to a two years old question of mine on MSE, I heard (for first time, I must confess) about two point Taylor series expansions; the commenter gave a link to a "quite" recent paper. Searching the Internet, I also found a rather old one.

However, I did not find papers showing examples on even basic functions. Does someone know where I could get illustrative examples ?

Being a very limited user of Mathematica (the coding of the method does not look simple - at least to me), I wonder if this method has been implemented in a standard package or as an add-on to Mathematica and if someone would accept to comment about his/her experience with it.

I hope that this question will not be closed too quickly for being out of topic.

Well, I'm not sure this is what you're looking for, but the first (older, Estes/Lancaster 1966) paper seems to be developing a way to calculate the Hermite interpolation at two points with the first k derivatives specified. The Hermite interpolation problem is solved by InterpolatingPolynomial, but not with the formulas in the paper:

func = Exp[Cos[1 + #]] &;
ip = Block[{k = 3, a = 0, b = 1, f = func}, (* general code for the problem *)
InterpolatingPolynomial[
{{{a}, Sequence @@ (NestList[D[#, x] &, f[x], k] /. x -> a)}, (* derivatives 0-k at a *)
{{b}, Sequence @@ (NestList[D[#, x] &, f[x], k] /. x -> b)}},(* derivatives 0-k at b *)
x]
];
{ip, ip /. x -> 1 + h} // Expand // N
{Series[func[x], {x, 0, 3}], Series[func[1 + h], {h, 0, 3}]} // Normal // N
(*
{1.71653 - 1.44441 x + 0.143992 x^2 + 0.460485 x^3 - 0.233256 x^4 - 0.0242814 x^5 + 0.0540343 x^6 - 0.013509 x^7,
0.659583 - 0.599757 h + 0.409921 h^2 - 0.107483 h^3 - 0.0169637 h^4 + 0.0162353 h^5 - 0.0405287 h^6 - 0.013509 h^7}

{1.71653 - 1.44441 x + 0.143992 x^2 + 0.460485 x^3,
0.659583 - 0.599757 h + 0.409921 h^2 - 0.107483 h^3}
*)


This code gives the polynomial interpolant, but possibly it was code for the coefficients that was desired. One way to get the coefficients:

SolveAlways[
ip == Sum[a[j] x (x (x - 1))^j, {j, 0, 3}] + Sum[b[j] (x - 1) (x (x - 1))^j, {j, 0, 3}],
{x}]


For the second (newer, López/Temme 2002) paper, below is an implementation of eqns. (3), (10) and (11). One can verify that the InterpolatingPolynomial code above produces the same polynomials as the code below; however, it will be in the form given in the paper. The function aLT[f, {z, z1, z2}, n] computes the coefficient $a_n(z_1, z_2)$ of López/Temme.

(* Lopez/Temme 2002, eqns. (10),(11) *)
ClearAll[aLT, aLTTerm, twoPtTaylor];
aLT[f_, {z_, z1_, z2_}, 0] := (f /. z -> z2)/(z2 - z1);
aLT[f_,                (* function expression in terms of symbol z *)
{z_, z1_, z2_},     (* symbol and two points *)
n_Integer?Positive  (* degree *)
] :=
Sum[((n + k - 1)! *
((-1)^(n + 1) n (D[f, {z, n - k}] /. z->z2) + (-1)^k k (D[f, {z, n - k}] /. z->z1))) /
(k! (n - k)! n! (z1 - z2)^(n + k + 1)),
{k, 0, n}];
(* Lopez/Temme 2002, eqn. (3) *)
aLTTerm[f_, {z_, z1_, z2_}, n_Integer?NonNegative] :=
(aLT[f, {z, z1, z2}, n] (z - z1) + aLT[f, {z, z2, z1}, n] (z - z2)) ((z - z1) (z - z2))^n;
twoPtTaylor[f_, {x_, x1_, x2_}, deg_] := Sum[aLTTerm[f, {x, x1, x2}, n], {n, 0, deg}]


Examples:

The coefficients:

Table[aLT[Sin[x], {x, 0, Pi}, k], {k, 0, 3}]
(*  {0, -(1/π^2), 1/π^4, -(2/π^6) + 1/(6 π^4)}  *)


Two examples from the MSE Q&A of the OP

twoPtTaylor[Sin[x], {x, 0, Pi}, 2]
(*  x^2 (-π + x)^2 (x/π^4 - (-π + x)/π^4) + x (-π + x) (-(x/π^2) + (-π + x)/π^2)  *)

Factor /@ twoPtTaylor[Sin[x], {x, -Pi, Pi}, 3]
(*
((π - x) x (π + x))/(2 π^2) +
(3 (π - x)^2 x (π + x)^2)/(8 π^4) -
((-15 + π^2) (π - x)^3 x (π + x)^3)/(48 π^6)
*)


The function aLTTerm gives a single term of eqn. (3) in López/Temme:

Grid[
Table[{k, aLTTerm[f[x], {x, x1, x2}, k]}, {k, 0, 2}],
Dividers -> All
]


The InterpolatingPolynomial approach will also work on symbolic input for the function f[x] and interval {x1, x2}; however, you need to pick a definite degree. One might try to obtain an expression for the coefficient for a general degree $n$ as a symbolic sum with something like this:

aLT[f0_, {z_, z1_, z2_}, n_ ] :=
With[{f = Evaluate[f0 /. z -> #] &},
Sum[((n + k - 1)! *
((-1)^(n + 1) n (Derivative[n - k][f][z2]) + (-1)^k k (Derivative[n - k][f][z1]))) /
(k! (n - k)! n! (z1 - z2)^(n + k + 1)),
{k, 0, n}]
];


Even with a specific function such Sin[x] or x^2, the Derivative[n - k][f] did not seem to simplify, and the Sum failed to return a closed form. (One can use SeriesCoefficient[f, {z, z1, n - k}] instead, but one has to supply some assumptions and judiciously apply simplification; moreover, it took much, much longer.) However, I did manage to get a result for $a_n(z_1, z_2)$ for $f(x)=e^x$ by manually injecting the derivative (although it's not what one traditionally calls a "closed form"):

a1 = aLT[Exp[x], {x, -1, 1}, n] /. Derivative[_][Function[__]] -> Exp


We can get the other coefficient $a_n(z_2, z_1)$ needed for the series too:

a2 = aLT[Exp[x], {x, 1, -1}, n] /. Derivative[_][Function[__]] -> Exp


We can then form the series approximation to $e^x$ and plot the error:

approx = Sum[(a1 (x + 1) + a2 (x - 1)) ((x + 1) (x - 1))^n, {n, 4}];
Plot[Exp[x] - approx, {x, -1, 1}]


• Thanks for the answer and clarification. But, then, this seems to mean that we should not get the formal expressions of the coefficients as in the comment to my post. Cheers. Apr 17, 2018 at 11:18
• @ClaudeLeibovici I think we can get formal expressions, at least for an explicit degree, for a generic f[x] and generic interval. Forgive me, I know the link has an anchor # to the comment, but when I click on it, MathJax makes the page dance around and does not land on the comment. I think it is the comment that links to the second (2002) paper. Do you mean you want code for formulas (10)-(11) in that paper? Apr 17, 2018 at 13:27
• What they give in the paper of 2002. I should be really amazed to see that working. Again, this is pure curiosity. I started with Mathematica version 2 or 3 long time ago and, being quite old, they have been decades without using it. May I confess that I am a fanatic lover of function approximation. So, if you are ready to waste some time, any example would be fantastic. Thanks in advance. Cheers. Apr 17, 2018 at 13:32
• @ClaudeLeibovici I added some more based on the 2002 paper. I'm still not sure just what sorts of examples you would like to produce. Apr 18, 2018 at 12:35
• This is very fine and I thank you very much for your work ! Cheers. Apr 19, 2018 at 4:42

Recently I wrote a recursive function to calculate the two-point Taylor expansion. You can find it here. The notebook includes a couple of examples that can get you started to play with it.

When I came up with this idea, I must admit that I didn't know the two-point Taylor expansion was a thing (I was playing with polynomial division). So, I didn't follow any paper in particular, which actually made my life simpler because I could write the expansion in what I consider was the simplest way: a recursion.

López/Temme's formula could be easily extended as laurent expansion and example consider the function:

-(1/2)+x/(2 \[Pi])-x/((-1+\[Pi]) \[Pi] (2+\[Pi]))+(\[Pi] x Cos[1])/(3 (-1+\[Pi]) (-1+x))-(x^2 Cos[1])/(3 (-1+\[Pi]) (-1+x))+(\[Pi] x Cos[2])/(6 (2+\[Pi]) (2+x))-(x^2 Cos[2])/(6 (2+\[Pi]) (2+x))

f[x_] := Cos[x]/((x + 2) (x - 1))


as López/Temme's formula $$-\frac{x^3}{4 \pi ^2}-\frac{x^3}{\pi ^3}+\frac{x^3}{(\pi -1)^2 \pi ^2 (2+\pi )}+\frac{2 x^3}{(\pi -1) \pi ^3 (2+\pi )}+\frac{x^3}{(\pi -1) \pi ^2 (2+\pi )^2}+\frac{x^2}{2 \pi }+\frac{3 x^2}{2 \pi ^2}-\frac{x^2}{(\pi -1)^2 \pi (2+\pi )}-\frac{3 x^2}{(\pi -1) \pi ^2 (2+\pi )}-\frac{x^2}{(\pi -1) \pi (2+\pi )^2}-\frac{x}{4}-\frac{1}{2}$$ as laurent series two point 0,PI $$-\frac{x^2 \cos (1)}{3 (\pi -1) (x-1)}-\frac{x^2 \cos (2)}{6 (2+\pi ) (x+2)}+\frac{x}{2 \pi }-\frac{x}{(\pi -1) \pi (2+\pi )}+\frac{\pi x \cos (1)}{3 (\pi -1) (x-1)}+\frac{\pi x \cos (2)}{6 (2+\pi ) (x+2)}-\frac{1}{2}$$

• An answer employing Mathematica code would be preferable. Aug 25, 2022 at 22:23
• Could you provide a code for $\sin(x)$ for $x \in (0,\pi)$ for any given number of terms. This could help m a lot. Thanks & cheers :-) Sep 9, 2022 at 9:21
• Hello, @Claude LeiboviciI did not understand your question very well. You have the expansion you need above using the formulas of Michael E2 Sep 11, 2022 at 7:45
• @Claude Leibovici If it helps you, I'll send you the expansion for cot for x (0 ,Pi) and from here you can get the formula you're looking for in product form. $$\sum _{k=2}^{\infty } \frac{(-1)^{-k} (\pi -x) x \cos (\pi k)}{\pi ^2 (k-1) k (\pi k-x)}-\sum _{k=2}^{\infty } \frac{(-1)^{-k} (\pi -x) x \cos (\pi k)}{\pi ^2 k (k+1) (\pi k+x)}-\frac{1}{2} x \left(\frac{4}{\pi ^2-x^2}+\frac{1}{\pi ^2}\right)+\frac{(\pi -x) (x+\pi )}{\pi ^2 x}=\cot (x)$$ Sep 11, 2022 at 20:04