Re: inverse command of Series[]
- To: mathgroup at smc.vnet.net
- Subject: [mg8081] Re: inverse command of Series[]
- From: Daniel Lichtblau <danl>
- Date: Tue, 5 Aug 1997 03:22:49 -0400
- Organization: Wolfram Research, Inc.
- Sender: owner-wri-mathgroup at wolfram.com
Matthias Weber wrote:
>
> In article <5s0v3n$q7n$1 at dragonfly.wolfram.com>, "Nguyen N. Anh"
> <anh at chm.ulaval.ca> wrote:
>
> > Hello,
> >
> > I've written a program in Mathematica and it works quite well. However, after
> > convergence it gave as result the following polinomial:
> >
> > 1 + x + x^2/2 + x^3/6 + x^4/24 + x^5/120 + x^6/720 + x^7/5040
> >
> > I know that this is the power series expansion of function E^x but I want that
> > my program recognizes it automatically and outputs the generating function of
> > the polinomial calculated ( E^x in this case).
> >
> > Is there any command or package which can do that ?
> >
> > Any help will be appreciated
> >
> > Regards
> >
>
> Here is a suggestion how to 'recognize' a polynomial as a
> truncated power series:
>
> Let us assume that the function you want to recognize
> satisfies some linear ode of order n. You can get its coefficients approximately
> by differentiating the polynomial n times and finding the coefficients
> of a linear combination which makes this zero.
> Then you can solve this ode explicitly, giving you the required function.
>
> For other types of functions, you might be able to modify this
> approach.
>
> However, you should be aware that this procedure will not
> prove anything.
>
> Good luck,
>
> Matthias Weber
I took a stab at this on the relatively simple example specified.
Possibly some details are not quite right, though theresult does work
out to what it should.
Say we suspect that the generating function satisfies an ode of order no
greater than four. We will take derivatives of order zero through four,
set up an equation in undetermined coefficients (the first we can
arbitrarily set to 1) for the powers of x, and solve it for the
indeterminates. Any indeterminate that is unsolved we arbitrarily set to
zero; this will give the simplest possible ode.
f[x_] = 1 + x + x^2/2 + x^3/6 + x^4/24 + x^5/120 + x^6/720 + x^7/5040;
(* Get derivates of orders 0 through 4 *)
derivs = Table[D[f[x],{x,j}], {j,0,4}]
(* Create array of undetermined coefficients; can fix one so set first
to 1 *)
coeffs = Array[a, 5, 0]; a[0] = 1;
eqn = Normal[coeffs . derivs + O[x]^4]
sol = SolveAlways[eqn == 0, x]
We get the solution {{a[1] -> -1 - a[2] - a[3] - a[4]}. Now set up our
ode. Plug in the values found for the a[j]'s. For any j>=1 for which we
found no value, use zero.
yderivs = Table[D[y[x], {x,j}], {j,0,4}]
deq = coeffs . yderivs /. sol[[1]] /. Thread[Rest[coeffs]->0]
Now solve the resulting ode.
DSolve[{deq==0, y[0]==1}, y[x], x]
yields {{y[x] -> E^x}}.
Here is a question. What should one do if the coefficients of the
approximating polynomial are only approximately correct? If I use a very
small perturbation, say
f[x_] = Table[1/(j!*(1+.000001*Random[Real,{-1,1}])), {j,0,7}] .
Table[x^j, {j,0,7}]
and then use
DSolve[Rationalize[{deq==0, y[0]==1, y[1]==f[1],
y[2]==f[2], y[3]==f[3]}], y[x], x]
I get a solution that, after massaging, looks a bit like E^x. But not
very much.
In[46]:= % /. Power[E,a_]:>Power[E,N[a]] // Chop // InputForm
Out[46]//InputForm=
{{y[x] -> 0.0714340061864687/E^(1.002577935612982*x) -
(0.02614998595252538 - 0.031089046185748034*I)*
E^((0.24497606817758646 - 0.9502216694997764*I)*x) -
(0.02614998595252552 + 0.031089046185747912*I)*
E^((0.24497606817758646 + 0.9502216694997764*I)*x) +
0.9808659657185822*E^(1.0000000284146515*x)}}
I'd be curious to see a better way.
Daniel Lichtblau
Wolfram Research
danl at wolfram.com