Plotting Taylor estimation with exact value

I'm trying to evaluate f(x) at x=-20 and plotting the Taylor estimation with the actual value. My results do not match with my hand calculations and I can't find the problem in my recursive for-loop.

%Using Taylor polynomials to evaluate f(x)=e^x at x=-20
n = 10;
x = -20;
%terms = 0;
terms = zeros(1, n+1);
for j = 0:1:n
    terms(j+1)=(x.^j)/factorial(j);
end;

%display(terms)
termsSum = sum(terms); %holds the estimate

display(terms)

display('The approximate estimate at x = -20')
display(termsSum)

%true value at x=-20
display('TRUE VALUE')
display(exp(x))

%absolute error between the approximated and true value
display('The error: ')

t = exp(x)
y = termsSum.^sign(x)
h = abs(t-y)
display(h)

%plot the approximation and exact graph
plot(-22:0.001:-18,termsSum,'o')
plot(x,exp(x),'*')

[John D'Errico]

Please learn to format your code so it is readable. Al that takes is to select it when you paste it in, then click on the "{} Code" button. Without that, pasted in code is an unreadable mess. I'll do it for you once.

[John D'Errico]

By the way, PLEASE DO NOT KEEP POSTING THE SAME QUESTION. If the first time you asked this, you did not get a clear enough answer, then prompt the person who answered it to clarify their response.

Answer 1

[John D'Errico]

As far as what is wrong, you have not said why you think there to be a problem.

When you write code like this, LOOK AT THE INTERMEDIATE RESULTS. THINK ABOUT WHAT IS HAPPENING.

Here is your terms vector. After 10 terms, we see a sequence of terms with alternating signs. This is a recipe for massive subtractive cancellation, which is surely what this assignment is designed to teach about.

terms =
          1          -20          200      -1333.3       6666.7       -26667        88889  -2.5397e+05   6.3492e+05  -1.4109e+06   2.8219e+06

I would expect this computation to have problems, at least if you want to get the correct answer. But the problems are entirely numerical. So lets look at what happens if I use a higher precision. I'll use essentially your code, but the computations will be done using my own HPF tool.

DefaultNumberOfDigits 100
terms = zeros(1,101,'hpf');
X = hpf(-20);
x = hpf(-20);
for i = 0:100
  terms(i+1) = x.^i./factorial(hpf(i));
end
>> terms
terms =
HPF array of size: 1  101
|1,1| 1
|1,2| -20
|1,3| 200
|1,4| -1333.333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333
|1,5| 6666.666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666667
|1,6| -26666.66666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666667
|1,7| 88888.88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888888889
|1,8| -253968.2539682539682539682539682539682539682539682539682539682539682539682539682539682539682539682540
|1,9| 634920.6349206349206349206349206349206349206349206349206349206349206349206349206349206349206349206349
|1,10| -1410934.744268077601410934744268077601410934744268077601410934744268077601410934744268077601410934744
|1,11| 2821869.488536155202821869488536155202821869488536155202821869488536155202821869488536155202821869489
|1,12| -5130671.797338464005130671797338464005130671797338464005130671797338464005130671797338464005130671797
|1,13| 8551119.662230773341884452995564106675217786328897440008551119662230773341884452995564106675217786329
|1,14| -13155568.71112426667982223537779093334648890204445760001315556871112426667982223537779093334648890204
|1,15| 18793669.58732038097117462196827276192355557434922514287593652673017752382831747911112990478069843149
|1,16| -25058226.11642717462823282929103034923140743246563352383458203564023669843775663881483987304093124199
|1,17| 31322782.64553396828529103661378793653925929058204190479322754455029587304719579851854984130116405249
|1,18| -36850332.52415760974740121954563286651677563597887282916850299358858338005552446884535275447195770881
|1,19| 40944813.91573067749711246616181429612975070664319203240944777065398153339502718760594750496884189868
|1,20| -43099804.12182176578643417490717294329447442804546529727310291647787529831055493432205000523035989334
|1,21| 43099804.12182176578643417490717294329447442804546529727310291647787529831055493432205000523035989334
|1,22| -41047432.49697311027279445229254566028045183623377647359343134902654790315290946125909524307653323175
|1,23| 37315847.72452100933890404753867787298222894203070588508493759002413445741173587387190476643321202887
|1,24| -32448563.23871392116426438916406771563672081915713555224777181741229083253194423814948240559409741641
|1,25| 27040469.36559493430355365763672309636393401596427962687314318117690902710995353179123533799508118034
|1,26| -21632375.49247594744284292610937847709114721277142370149851454494152722168796282543298827039606494427
|1,27| 16640288.84036611341757148162259882853165170213186438576808811149348247822150986571768328492004995713
|1,28| -12326139.88175267660560850490562876187529755713471435982821341592109813201593323386495058142225922751
|1,29| 8804385.629823340432577503504020544196641112239081685591581011372212951439952309903536129587328019646
|1,30| -6071990.089533338229363795520014168411476629130401162476952421636008932027553317174852503163674496308
|1,31| 4047993.393022225486242530346676112274317752753600774984634947757339288018368878116568335442449664205
|1,32| -2611608.640659500313704858288178136951172743712000499990087063069251153560237985881656990608032041423
|1,33| 1632255.400412187696065536430111335594482964820000312493804414418281970975148741176035619130020025889
|1,34| -989245.6972195076945851735940068700572624029212123106023057057080496793788780249551731025030424399329
|1,35| 581909.2336585339379912785847099235630955311301248885895915915929703996346341323265724132370837881958
|1,36| -332519.5620905908217093020484056706074831606457856506226237666245545140769337899008985218497621646833
|1,37| 184733.0900503282342829455824475947819350892476586947903465370136413967094076610560547343609789803796
|1,38| -99855.72435152877528807869321491609834329148522091610289002000737372795103116813840796451944809750250
|1,39| 52555.64439554146067793615432364005175962709748469268573158947756511997422693059916208658918320921184
|1,40| -26951.61251053408239894161760186669321006517819727830037517409105903588421893876880107004573497908300
|1,41| 13475.80625526704119947080880093334660503258909863915018758704552951794210946938440053502286748954150
|1,42| -6573.564026959532292424784780943095904893945901775195213457095380252654687546041170992694081702215365
|1,43| 3130.268584266443948773707038544331383282831381797712006408140657263168898831448176663187657953435888
|1,44| -1455.938876402997185476142808625270410829223898510563723910763096401473906433231710075901236257412041
|1,45| 661.7903983649987206709740039205774594678290447775289654139832256370335938332871409435914710260963822
|1,46| -294.1290659399994314093217795202566486523684643455684290728814336164593750370165070860406537893761699
|1,47| 127.8822025826084484388355563131550646314645497154645343795136667897649456682680465591481103432070304
|1,48| -54.41795854579082912290874736730002750275087221934661037426113480415529602905023257836089801838597038
|1,49| 22.67414939407951213454531140304167812614619675806108765594213950173137334543759690765037417432748766
|1,50| -9.254754854726331483487882205323133929039263982882076594262097755808723814464325268428724152786729656
|1,51| 3.701901941890532593395152882129253571615705593152830637704839102323489525785730107371489661114691862
|1,52| -1.451726251721777487605942306717354341810080624765815936354838863656270402268913767596662612201839946
|1,53| 0.5583562506622221106176701179682132083884925479868522832133995629447193854880437567679471585391692100
|1,54| -0.2107004719480083436293094784785710220333934143346612389484526652621582586747334931199800598261015887
|1,55| 0.07803721183259568282567017721428556371607163493876342183276024639339194765730870115554817030596355136
|1,56| -0.02837716793912570284569824625974020498766241270500488066645827141577889005720316405656297102035038231
|1,57| 0.01013470283540203673060651652133578749559371882321602880944938264849246073471541573448677536441085083
|1,58| -0.003556036082597205870388251410995013156348673271303869757701537771400863415689619555960272057688017834
|1,59| 0.001226219338826622713926983245170694191844370093553058537138461300483056350237799846882852433685523391
|1,60| -0.0004156675724836009199752485576849810819811424045942571312333767120281546949958643548755431978594994546
|1,61| 0.0001385558574945336399917495192283270273270474681980857104111255706760515649986214516251810659531664849
|1,62| -0.00004542814999820775081696705548469738600886802236002810177413953136919723442577752512301018555841524094
|1,63| 0.00001465424193490572606998937273699915677705420076129938766907726818361201110508952423322909211561781966
|1,64| -0.000004652140296795468593647419916507668818112444686126789736215005772575241620663341026421934004958037986
|1,65| 0.000001453793842748583935514818723908646505660138964414621792567189303929763006457294070756854376549386871
|1,66| -0.0000004473211823841796724660980688949681555877350659737297823283659396706963096791674063867244235536574987
|1,67| 0.0000001355518734497514158988175966348388350265863836284029643419290726274837302058083049656740677435325754
|1,68| -0.00000004046324580589594504442316317457875672435414436668745204236688735148768065845024028826091574433808220
|1,69| 0.00000001190095464879292501306563622781728138951592481373160354187261392690814137013242361419438698362884770
|1,70| -0.000000003449552072113891308134967022555733736091572409777276388948583746929896049313745975128807821341694987
|1,71| 0.0000000009855863063182546594671334350159239245975921170792218254138810705513988712324988500368022346690557105
|1,72| -0.0000000002776299454417618759062347704270208238303076386138653029334876255074363017556334788836062632870579466
|1,73| 0.00000000007711942928937829886284299178528356217508545517051813970374656264095452826545374413433507313529387406
|1,74| -0.00000000002112861076421323256516246350281741429454396032068990128869768839478206253848047784502330770829969152
|1,75| 0.000000000005710435341679252044638503649410111971498367654240513861810186052643800686075804822979272353594511223
|1,76| -0.000000000001522782757781133878570267639842696525732898041130803696482716280705013516286881286127805960958536326
|1,77| 0.0000000000004007323046792457575184914841691306646665521160870536043375569159750035569176003384546857791996148227
|1,78| -0.0000000000001040863129037001967580497361478261466666369132693645725552095885649359888097663216765417608310687851
|1,79| 0.00000000000002668879818043594788667941952508362735041972135111912116800245860639384328455546709654916944386379105
|1,80| -0.000000000000006756657767198974148526435322805981607701195278764334472912014837061732477102649897860549226294630646
|1,81| 0.000000000000001689164441799743537131608830701495401925298819691083618228003709265433119275662474465137306573657662
|1,82| -0.0000000000000004170776399505539597855824273337025683766169925163169427723465948803538566112746850531203226107796695
|1,83| 0.0000000000000001017262536464765755574591286179762361894187786625163275054503889952082577100669963544195908806779682
|1,84| -0.00000000000000002451235027625941579697810328144005691311295871385935602540973228800198980965469791672761226040432968
|1,85| 0.000000000000000005836273875299860904042405543200013550741180646156989529859460068571902335632070932554193395334364210
|1,86| -0.000000000000000001373240911835261389186448363105885541350866034389879889378696486722800549560487278248045504784556285
|1,87| 0.0000000000000000003193583515895956719038252007222989631048525661371813696229526713308838487349970414530338383219898336
|1,88| -0.00000000000000000007341571300910245331122418407409171565628794623843249876389716582319168936436713596621467547631950199
|1,89| 0.00000000000000000001668538932025055757073276910774811719461089687237102244634026495981629303735616726504878988098170500
|1,90| -0.000000000000000000003749525689943945521512981821965869032496830757836184819402306732542987199405880284280626939546450561
|1,91| 0.0000000000000000000008332279310986545603362181826590820072215179461858188487560681627873304887568622853956948754547667913
|1,92| -0.0000000000000000000001831270178238801231508171830019960455431907574034766700562787170961165909355741286583944781219267673
|1,93| 0.00000000000000000000003981022126606089633713417021782522729199799073988623262093015589046012846425524536052053872215799290
|1,94| -0.000000000000000000000008561337906679762653147133380177468234838277578470157552888205567840887841775321582907642735947955462
|1,95| 0.000000000000000000000001821561256740375032584496463867546432944314378397905862316639482519337838675600336788860156584671375
|1,96| -0.0000000000000000000000003834865803663947437019992555510624069356451322942959710140293647409132291948632287976547698072992368
|1,97| 0.00000000000000000000000007989303757633223827124984490647133477825940256131166062792278432102358941559650599951141037652067434
|1,98| -0.00000000000000000000000001647279125285200789097934946525182160376482527037353827379851223113888441558690845350750729412797409
|1,99| 3.361794133235103651220275401071800327298943932729293525265002496150792737874879276226021896760811039e-27
|1,100| -6.791503299464855861051061416306667327876654409554128333868691911415742904797735911567721003557194018e-28
|1,101| 1.358300659892971172210212283261333465575330881910825666773738382283148580959547182313544200711438804e-28

See that the first few terms I computed are the same as what you had, which is what one should expect.

sum(terms)
ans =
0.000000002061153622438557827988422171469799527665652835635718872481909232440581903875866301133659586488599105

exp(x)
ans =
0.000000002061153622438557827965940380155820976375807275599103692972244661629164023784559353279910927905581367

The problem is that the individual terms are massive. They don't drop below 1 in magnitude until about the 53rd term. At the end, they are pretty small, so the total sum is reasonable. The two agree out to about 29 decimal digits, as we should expect by looking at the individual terms.

sum(terms) - exp(x)
ans =
2.248179131397855128984556003661517950966457081141788009130694785374865858301773770841002000000000000e-29

So, there is no real problem in what you did, except that you did not go far enough. Be careful though. IF you try to use double precision to compute things like the factorial of a very large number, at some point it will overflow into an inf.

[John D'Errico]

It is important to recognize that just because the final term in that sum is on the order of 1e-28, that the result need not be truly be accurate to roughly that many digits. So stopping as soon as a term drops below some tolerance can still be somewhat dangerous.

For example, consider the alternating sum of the inverse of the integers, 1:N.

format long g
N = 1:1000;s=mod(N,2)*2 - 1;sum(1./N.*s)
ans =
         0.692647430559822

N = 1:10000;s=mod(N,2)*2 - 1;sum(1./N.*s)
ans =
         0.693097183059958

N = 1:100000;s=mod(N,2)*2 - 1;sum(1./N.*s)
ans =
          0.69314218058496

N = 1:1000000;s=mod(N,2)*2 - 1;sum(1./N.*s)
ans =
          0.69314668056024

N = 1:10000000;s=mod(N,2)*2 - 1;sum(1./N.*s)
ans =
         0.693147130560056

N = 1:100000000;s=mod(N,2)*2 - 1;sum(1./N.*s)
ans =
         0.693147175560058

The final term in each of those sums is pretty small, but as you see, the series is very slowly convergent. The final value, when this (alternating harmonic) series does converge, is actually

log(2)
ans =
       0.693147180559945

As you can see, the problem here is not so much massive subtractive cancellation, but just a very slowly convergent series. Even after 100 million terms, the sum is still pretty far away.

[cee878]

I see what you mean. I was also trying to look at it graphically as well but I think there's syntax error in my plot of the 2 measurements I want. However, I'm not sure what's wrong with it.

[John D'Errico]

Looking at the latter part of your script ... so then what was the purpose of this line?

y = termsSum.^sign(x)

x was -20, as used in your series solution. That will yield a reasonable value as I showed, if the series was run far enough, and you used a sufficiently high precision. sign(-20) will be -1, so that will invert what was your answer. Perhaps this is the problem of which you speak.