7
$\begingroup$

I have some complicated functions from numerical procedures that go to simple trig functions in the limit of small time-steps. For example, I can calculate that one of these functions goes to:

x - 2x^3 / 3 + 2x^5 / 15

This is $\sin{x} \cos{x}$ minus terms of $O(x^7)$ and higher. Is there any way for Mathematica to recognize such a pattern automatically? I have other similar series expansions and I don't want to have to find all of them manually.

---Edit---

Thank you Szabolcs and corey979! You asked for more examples/background so here it is. I wanted to convert a 2x2 matrix of polynomials to trigonomic functions. I can obtain this matrix to arbitrary order in some small parameter h. To first order that looks like:

{{1 + h (t - (2 t^3)/3 + (2 t^5)/15), h (t^2 - t^4/3 + (2 t^6)/45)}, {h (t^2 - t^4/3 + (2 t^6)/45), 1 + h (-t + (2 t^3)/3 - (2 t^5)/15)}}

To second order that looks like:

{{1 + h^2 (t^2/2 - t^4/6) + h (t - (2 t^3)/3 + (2 t^5)/15), h (t^2 - t^4/3 + (2 t^6)/45)}, {h (t^2 - t^4/3 + (2 t^6)/45), 1 + h^2 (t^2/2 - t^4/6) + h (-t + (2 t^3)/3 - (2 t^5)/15)}}

etc

For a single matrix element, I was able to collect the coefficients of various orders of h and use With the data and FindFunction approach to turn those coefficents into trigonomic functions. It seems to work okay if I have the coefficients have enough terms, but I can get pretty wacky numbers if it doesn't have enough terms. I use this code:

findtrigfunc[inseries_] :=  Module[{outseries=inseries},
outseries = outseries //CoefficientList[#,h]&;
For[i=1, i <= Length[outseries],i++, 
{data = Table[Evaluate@{x, outseries[[i]] //. t->x}, {x, 0, 1, 0.01}];
outseries[[i]] = FindFormula[data, x, TargetFunctions -> {Times,Sin, Cos}]}];
Sum[outseries[[i]]*h^(i- 1) , {i,1,Length[outseries]}]//.x->t]

Here's what that produces for the first matrix element up to second order in h and order 7 in t:

1. + 1.0043 h Cos[t] Sin[t] + h^2 (-0.00618535 + 0.596221 Sin[Sin[t]]^2)

When I increase the order in t from 7 to 11, I get extra terms that get me a more reasonable answer:

1. + h^2 (0.49996 - 0.499954 Cos[t]^2) + 1.00001 h Cos[t] Sin[t]

$\endgroup$
  • $\begingroup$ If the expressions appear exactly as written, you should be able to use ReplaceAll (i.e. /. replacement rules); it seems difficult to do in more general cases, though. $\endgroup$ – MarcoB Dec 7 '16 at 19:19
  • $\begingroup$ Can you provide a few more examples of such series, and what functions do they come from? $\endgroup$ – corey979 Dec 7 '16 at 19:31
6
$\begingroup$

If you have sufficiently many terms, FindGeneratingFunction may be able to do exactly this. But I wasn't having much luck with it, so I went in a roundabout way:

This is the expression. The more terms you have the better it works.

expr = Normal@Series[Sin[x] Cos[x], {x, 0, 8}]
(* x - (2 x^3)/3 + (2 x^5)/15 - (4 x^7)/315 *)

Extract coefficients:

coeff = CoefficientList[expr, x]
(* {0, 1, 0, -(2/3), 0, 2/15, 0, -(4/315)} *)

Try FindGeneratingFunction. It doesn't work here. But if we had computed the series to order 12, it would work, and it would return 1/2 Sin[2 x] (try it).

FindGeneratingFunction[coeff, x]
(* FindGeneratingFunction[{0, 1, 0, -(2/3), 0, 2/15, 
  0, -(4/315)}, x] *)

Construct sequence and get rid of the factorials, in the hope that it leads to simpler terms:

seq = {#1, #1! #2} & @@@ Transpose@{Range[0, Length[coeff] - 1], coeff}
(* {{0, 0}, {1, 1}, {2, 0}, {3, -4}, {4, 0}, {5, 16}, {6, 0}, {7, -64}} *)

Now try FindSequenceFunction. When this method fails, this is where it will fail.

fun = FindSequenceFunction[seq]
(* -2^(-2 + #1) ((-I)^(1 + #1) + I^(1 + #1)) & *)

Re-sum the result:

Sum[fun[k]/k! x^k, {k, 0, Infinity}]
(* -(1/4) I E^(-2 I x) (-1 + E^(4 I x)) *)

Simplify:

FullSimplify[%]
(* Cos[x] Sin[x] *)
$\endgroup$
  • $\begingroup$ If you consider that the length of the pattern in the series of Sin[] and Cos[] is 4, a series of order 8 is barely long enough to guess there is a pattern (at least naively). $\endgroup$ – Michael E2 Dec 7 '16 at 19:47
  • $\begingroup$ @MichaelE2 I did edit the post to note that it works with 12 (a few minutes ago). $\endgroup$ – Szabolcs Dec 7 '16 at 19:51
  • $\begingroup$ @MichaelE2 Do you know how FindGeneratingFunction works, or how such a thing may be implemented? $\endgroup$ – Szabolcs Dec 7 '16 at 21:10
  • 1
    $\begingroup$ No, I do not for sure, but I'd bet a small amount of money the magic word is "holonomic." For an overview, see Kauers, The Holonomic Toolkit, which I've only skimmed. Sect. 4 describes "guessing" for this sort of problem. $\endgroup$ – Michael E2 Dec 7 '16 at 21:37
5
$\begingroup$

A similar idea to Szabolcs, but I was waiting for a few more examples from the OP. Here goes.

f[x_] := x - 2 x^3/3 + 2 x^5/15

Create some data from it

data = Table[Evaluate@{x, Normal@f[x]}, {x, 0, 0.1, 0.01}];

and try to FindFormula:

FindFormula[data, x, TargetFunctions -> {Times, Sin, Cos}]

Cos[x] Sin[x]


A few made-up examples:

g[x_] := Normal@Series[Tan[x] + Cos[x], {x, 0, 5}]
g[x]

1 + x - x^2/2 + x^3/3 + x^4/24 + (2 x^5)/15

data2 = Table[Evaluate@{x, Normal@g[x]}, {x, 0, 0.1, 0.01}];

FindFormula[data2, x, TargetFunctions -> {Times, Plus, Sin, Cos, Tan}]
  1. Cos[x] + 1. Tan[x]

Also works fine.


h[x_] := Normal@Series[Sin[x] Tan[x], {x, 0, 5}]
h[x]

x^2 + x^4/6

data3 = Table[Evaluate@{x, Normal@h[x]}, {x, 0, 0.1, 0.01}];

FindFormula[data3, x, TargetFunctions -> {Times, Sin, Cos, Tan, Sec}]

0.999993 Sin[x] Tan[x]

Almost perfect (doesn't work without Sec though).


k[x_] := Normal@Series[Cos[x]^2 + Sin[x], {x, 0, 4}]
k[x]

1 + x - x^2 - x^3/6 + x^4/3

data4 = Table[Evaluate@{x, Normal@k[x]}, {x, 0, 0.1, 0.01}];

FindFormula[data4, x, TargetFunctions -> {Times, Plus, Sin, Cos}]

Cos[x]^2 + Sin[x]


Unfortunately, this fails

p[x_] := Normal@Series[Sin[x]^2 + Cos[x], {x, 0, 6}]
p[x]

1 + x^2/2 - (7 x^4)/24 + (31 x^6)/720

data5 = Table[Evaluate@{x, Normal@p[x]}, {x, 0, 0.1, 0.01}];

FindFormula[data5, x, TargetFunctions -> {Times, Plus, Sin, Cos}]

1 + 0.499998 x^2 - 0.29104 x^4

probably because Sin[x]^2 is small compared to Cos[x].

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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