Numerieke berekening van de afgeleide van de sinus

Remove["Global`*"]

(a) Schrijf een functie d1sin met argumenten x en h die een numerieke benadering berekent van de afgeleide van sin(x) met behulp van een eerste orde differentie quotiënt, dat wil zeggen f'(x)=(f(x+h)-f(x))/h.

d1sin[x_,h_]:=(Sin[x+h]-Sin[x])/h

(b) Schrijf nu een functie diff1sin met argumenten x en h die het verschil berekent tussen d1sin en de echte afgeleide (gebruik daarvoor een accent ' ) van de sinus, in het punt x, met differentie h.

diff1sin[x_,h_]:=d1sin[x,h]-Sin'[x]

(c) Maak een lijst hlist met tweemachten van -1 tot -64 (dwz 2^-1, 2^-2, ..., 2^-64)

hlist=Table[2^(-i),{i,1,64}]

(d) Pas de functie diff1sin toe op hlist in het punt x=1., vergeet niet de punt achter de 1.  Deze punt zorgt ervoor dat het double precision getal 1 wordt gebruikt, en niet de exacte integer 1. Noem de output o1.

o1=diff1sin[1.,hlist]

(e) Plot abs(o1) op logarithmische schaal mbv ListLogPlot. Gebruik de opties AxesLabel (bedenk zelf de juiste labels) en GridLines->Automatic. Noem het plaatje plote

plote=ListLogPlot[Abs[o1],AxesLabel->{"-log(2,h)","error"},GridLines->Automatic]

Bovenstaand plaatje heeft als voordeel dat de as index verwijst naar de bits van de mantisse, maar als nadeel dat de waarde van h niet direct valt af te lezen, wat het schatten van de hellingen bemoeilijkt. Daarom is een dubbellogplaatje makkelijker af te lezen. Maak daarvoor eerst een dubbele lijst met coordinaten door de functie Transpose te laten werken op een dubbele lijst van hlist en abs(o1). Gebruik dan de functie ListLogLogPlot met dezelfde opties als boven. Noem het plaatje plote2

plote2=ListLogLogPlot[Transpose[{hlist,Abs[o1]}],AxesLabel->{"h","error"},GridLines->Automatic]

(f) Schrijf een functie d2sin met argumenten x en h die een numerieke benadering berekent van de afgeleide van sin(x) met behulp van een tweede orde differentie quotiënt, dat wil zeggen f'(x)=(f(x+h)-f(x-h))/(2h).

d2sin[x_,h_]:=(Sin[x+h]-Sin[x-h])/(2h)

(g) Schrijf nu een functie diff2sin met argumenten x en h die het verschil berekent tussen d2sin en de echte afgeleide van de sinus, in het punt x, met differentie h.

diff2sin[x_,h_]:=d2sin[x,h]-Sin'[x]

(h) Pas de functie diff2sin toe op hlist in het punt x=1, vergeet niet de punt achter de 1. Noem de output o2

o2=diff2sin[1.,hlist]

(i) Plot abs(o2) op logarithmische schaal mbv ListLogPlot met een andere kleur. Noem het plaatje ploti

ploti=ListLogPlot[Abs[o2],PlotStyle->Red,GridLines->Automatic]

ploti2=ListLogLogPlot[Transpose[{hlist,Abs[o2]}],AxesLabel->{"h","error"},PlotStyle->Red,GridLines->Automatic]

(j) Laat de beide plaatjes samen zien (hint: gebruik Show[ploti,plote])

Show[ploti,plote]

Show[ploti2,plote2]

(k) Beschrijf en verklaar wat je ziet in de bovenstaande figuren (ongeveer 10 zinnen). Welke waarde van differentie h is optimaal voor de verschillende benaderingen?

(l) Wanneer de exacte integer 1 wordt gebruikt, dus zonder punt erachter, zal Mathematica proberen exact te rekenen. Daarbij komen een aantal waarschuwingen dat de $MaxExtraPrecision limiet bereikt wordt.
Becommentarieer het nieuwe plaatje. Combineer ook dit plaatje met de beide plaatjes uit onderdeel j.

plotl=ListLogPlot[Abs[diff2sin[1,hlist]],PlotStyle->Blue]

N[1,50]

plotl=ListLogPlot[Abs[diff2sin[N[1,50],hlist]],PlotStyle->Blue]

plotl=ListLogPlot[Abs[N[diff2sin[1,hlist],1]],PlotStyle->Blue]

(m) Schrijf nu een functie met argument x die Module gebruikt en waarin het combinatieplaatje van onderdeel l getekend wordt. Test de functie, en gebruik daarna Manipulate om een aantal x waarden tussen 0 en 2π te onderzoeken. Om te zorgen dat het argument exact is gebruik je bv {x,0,2π,π/128}.

m[y_,color_:Black]:=Module[{hlist,o1},d1sin[x_,h_]:=(Sin[x+h]-Sin[x])/h;diff1sin[x_,h_]:=d1sin[x,h]-Sin'[x];hlist=Table[2^(-i),{i,1,64}];o1=diff1sin[y,hlist];plote=ListLogPlot[Abs[o1],PlotStyle->color,AxesLabel->{"-log(2,h)","error"},GridLines->Automatic]]

Show[m[1,Red],m[1.]]

Becommentarieer het resultaat.

(n) Gebruik nu Manipulate om de volgende x waarden te onderzoeken x=1.+2kπ, met {k,0,10,1}.

Manipulate[m[1.+2k π],{k,0,10,1}]

Wat merk je op ? Is de precisie relatief of absoluut?

(o) Hoe beoordeel je de geschiktheid van Mathematica als PSE voor dit probleem?