MathGroup Archive 1997

[Date Index] [Thread Index] [Author Index]

Search the Archive

Re: inverse command of Series[]

  • To: mathgroup at smc.vnet.net
  • Subject: [mg8097] Re: inverse command of Series[]
  • From: Daniel Lichtblau <danl>
  • Date: Tue, 12 Aug 1997 00:54:38 -0400
  • Organization: Wolfram Research, Inc.
  • Sender: owner-wri-mathgroup at wolfram.com

The basic set-up is copied from my previous reply:

>f[x_] = 1 + x + x^2/2 + x^3/6 + x^4/24 + x^5/120 + x^6/720 + x^7/5040;
>derivs = Table[D[f[x],{x,j}], {j,0,4}]
>coeffs = Array[a, 5, 0]; a[0] = 1;
>eqn = Normal[coeffs . derivs + O[x]^4]
>sol = SolveAlways[eqn == 0, x]
>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.
> ...
>I'd be curious to see a better way.


I have seen the error of my ways. The polynomial approximation to E^x is
very good for the first several derivatives at the origin. It is not
good for values of f (or derivatives) away from the origin. This is
obvious, because the polynomial is (approximately) a truncation of the
power series, not a result of interpolation. So here is a better way.

dsol = DSolve[Rationalize[{deq==0, y[0]==f[0], y'[0]==f'[0],
        y''[0]==f''[0], y'''[0]==f'''[0]}], y[x], x];

ndsol = dsol /. Power[E,a_] :> Power[E,N[a]];

ndsol2 = ndsol /. Times[a_, Exp[b_]] :> N[a]*Exp[b] // Chop[#, 10^(-6)]&

I get

In[24]:= InputForm[ndsol2]
Out[24]//InputForm= {{y[x] ->
0.9999997405748506*E^(0.99999998089691*x)}}

I chose to Chop at 10(-6) because 10(-7) would not suffice to get rid of
all other terms. I should also mention that I am using the development
version of Mathematica. The versions on the market may not quite do
this. One must sometimes use Rationalize[expression,0] (explicitly give
tolerance of zero) because default behavior is a bit different.


Daniel Lichtblau
Wolfam Research
danl at wolfram.com


  • Prev by Date: Math 3.0 on AIX
  • Next by Date: Sorting Eigenvalues/EigenVectors
  • Previous by thread: Re: inverse command of Series[]
  • Next by thread: Combining ListPlots