Path: news.mathworks.com!newsfeed-00.mathworks.com!solaris.cc.vt.edu!news.vt.edu!news.glorb.com!news2.glorb.com!Xl.tags.giganews.com!border1.nntp.dca.giganews.com!nntp.giganews.com!local01.nntp.dca.giganews.com!nntp.supernews.com!news.supernews.com.POSTED!not-for-mail
NNTP-Posting-Date: Wed, 18 Mar 2009 18:13:31 -0500
Subject: Re: How to take 200th order derivative, and then evaluate it numerically?
Date: Wed, 18 Mar 2009 16:13:30 -0700
From: Ronald Bruck <bruck@math.usc.edu>
Newsgroups: comp.soft-sys.math.maple,sci.math.num-analysis,comp.soft-sys.matlab
Message-ID: <180320091613300494%bruck@math.usc.edu>
References: <f6ebcf2b-4782-4984-aebc-7cdb5b28ccee@p2g2000prf.googlegroups.com> <rbisrael.20090318173245$27f5@news.acm.uiuc.edu> <gprk0t$g8$1@aioe.org>
MIME-Version: 1.0
Content-Type: text/plain; charset=ISO-8859-1
Content-transfer-encoding: 8bit
User-Agent: Thoth/1.8.4 (Carbon/OS X)
Lines: 74
X-Trace: sv3-8TAwEOxLpWJE/qKEiXkk/97Wm+JwSNVQ66ek4zfM2m2hyZwjG28PId6At+JPt4lMN0mTjGASqgLUsgA!Xuf2oR2OmCsWndDhbS95gUU/nZ0cJljnU1XVjgole+a4UQhxj0rViFoFwbi8D27r/tr+yRO/rYg=
X-Complaints-To: www.supernews.com/docs/abuse.html
X-DMCA-Complaints-To: www.supernews.com/docs/dmca.html
X-Abuse-and-DMCA-Info: Please be sure to forward a copy of ALL headers
X-Abuse-and-DMCA-Info: Otherwise we will be unable to process your complaint properly
X-Postfilter: 1.3.39
Xref: news.mathworks.com comp.soft-sys.math.maple:24558 sci.math.num-analysis:107044 comp.soft-sys.matlab:526004


In article <gprk0t$g8$1@aioe.org>, dpb <none@non.net> wrote:

> Robert Israel wrote:
> ...
> > Perhaps what you might try is computing a Taylor series up to order x^200.  
> > There still might be some numerical difficulties, ... [snip for effect]
> 
> Think???  :)
> 
> >  For example, in Maple:
> > 
> >   f:= x^3*sin(x)/(x - cos(x));
> >   Digits:= 20;
> >   200! * coeff( series(f, x=0.3, 201), x-0.3, 200);
> > 
> >     -.90246797086526191165e446
> 
> Hmmmm....any idea how many of those digits are actually meaningful in 
> front of that 10**446 exponent???

I don't get ANY.  Using Mathematica, and the commands

f[x_] = x^3 * Sin[x]/(x - Cos[x]); 
a = 3/10;
q = 200! * Coefficient[Series[f[x + a], {x, N[a, 2000], 201}], x^200]

(i.e., using 2000-digit precision), I get as output

8.79706766309285531402058935431446364695557682377643868580539624307543\
3813635170485359021030321735306371237255460885001870669592062272255754\
6811798739911826558405041319246793190060013582632484660798804338639829\
3677497582276781551805425896909210803391220133566681155945000972282627\
0692882515065842765398727066878423388527454848431676355424447555976171\
7362630436351635192836439246091937861477700209751375257425678609745663\
1485232770451072262442515742371840311536121150067032833905408799375142\
8105502883510239941596597308186610428729501136776398590701695574041338\
6975242778619445713968545975844143146964882124308981769154558016449941\
1567681808863587798175197072693803457483601743383577315698634205703705\
3868859280352521965983155519440299399272535965564387743166865563501490\
1480369166792083169418839803175212971169119214904804916044561800575025\
0771513386032793328831890682203270760620393851840263530296059615280085\
7033088303440083093200749796467206434398433793567452930766246051370380\
3698107871076710080131452293836149387506048965829763565001806694499146\
1810219474350028345839056832965424219228595190047972663418201733913850\
8246700095809447655044553220714537798969707319935194914526012068081080\
6210154194150405075121029247671378458058680688604742558715468066043142\
2741282672261152073353657296454567641371140261797813339278678807753088\
0168735036976654246639887730891366191064192337684932419310928233540923\
9981278936062880224390643444860496474151019886455346989575020619857174\
9695402221976286543825694523159480697218429325203913636735567724996871\
2453332488236538342660984769811105037613490027917909424973805659192699\
3628345493851270616604909474801236609765889934117075983314385416221233\
9858564165392993939159640504400071083456475018607584848375789700395435\
3978324345557471056134749552317765853161430754113665525334594231640545\
9493374689336586922233202178132503930023983656978124892451476942480899\
7644028138163525412702361745808942154119575168496464593462586541700193\
5315806189348358721672465614921305949*10^548

(that has 1996 digits, I believe).  Which isn't very far from what
Mathematica returns with 20-digit precision,

8.797067663092855*10^548

The exact answer is a rational function of cos(3/10) and sin(3/10). 
The easiest way to find the coefficients might be to set up a
recurrence formula for them.  Ugly.  But don't ask Mathematica to
return the exact answer; it will take forever, for a simple demo
example.

As for the OP, I would suggest doing this with 16 digits, then 32, then
64, ... and compare the answers.  Crude, but it usually works.

-- 
Ron Bruck