Kvadratur

Kvadratur och interpolation

Matlabs integral är en bra rutin, men den har svårt med ändpunkts-singulariteter.

>> integral(@(x) 1 ./ x.^0.7, 0, 1)
Warning: Infinite or Not-a-Number value encountered.
> In funfun/private/integralCalc>iterateScalarValued at 350
In funfun/private/integralCalc>vadapt at 133
In funfun/private/integralCalc at 76
In integral at 88
ans =
Inf

Man tvinga fram ett svar genom att öka den absoluta toleransen. Jag har subtraherat det exakta svaret, 1/0.3.

>> integral(@(x) 1 ./ x.^0.7, 0, 1, 'AbsTol', 1e-4) - 1 / 0.3
ans =
-3.0286e-05

För att klara av svåra fall kan man använda variabelsubstitution eller partiell integration för att bli av med singulariteten och få en snäll integrand som integral kan hantera. I följande uppgift får du träna på denna teknik. Integrationsintervallet är [0, 1] i samtliga tre fall.

a) Använd partiell integration följd av integral för att approximera den bestämda integralen av cos(x) / x^0.9.
b) Använd variabelsubstitution följd av integral för att approximera den bestämda integralen av cos(x) / x^0.9. Ledning: testa ansatsen t = x^p där p är ett reellt tal.

c) Frivillig: Använd variabelsubstitution följd av integral för att approximera den bestämda integralen av log(x) / (x^0.85 + x^0.9 + x^0.95). Ledning: testa ansatsen t = x^p där p är ett reellt tal.

Länkarna i det första stycket behöver du bara läsa om intresse föreligger.
I följande uppgift kommer du att ta fram formen på ett bord. Det är en typ av problem som jag har haft flera gånger. Under de senaste åren har vi tagit fram ett bord som haft formen av en superellips, ett mellanting mellan en ellips och en rektangel. Formen, som har använts för fontäner och andra arkitektoniska skapelser, var ett påfund av den danske mångsysslaren Piet Hein, en länk till (en del av en utmärkt matematisk uppslagsbok), slutligen Wikipedias artikel. Här kan du se ett bord, Supercirkel av Bruno Mathsson & Piet Hein.

Eftersom som jag har lagt till ovanstående två övningar i år har jag förenklat designen av bordet. I år kommer du att bestämma utseendet på ett bord som har formen av en epicykloid. Vi kommer att generalisera bilderna från länken till vänster, man talar om shortened eller lengthened epicykloid (kanske förkortad, förlängd på svenska). Här är tre epicykloider (de tre svarta kurvorna).

Den blå cirkeln rullar utmed den röda, utan att slira. Den svarta ringens väg genererar kurvan. I nästa bild har den svarta punkten flyttats in utmed den blå cirkelns radie, r. Ringen befinner sig på avståndet 0.5 r från den blå cirkelns centrum (shortened). I den tredje bilden har ringen flyttats ut och befinner sig på avståndet 2 r från blå cirkelns centrum (lengthened). Låt q vara graden av förkortning eller förlängning, så i bilderna är q = 1, q = 0.5 respektive q = 2. De tre bilderna är inte ritade i samma skala.

Om R är radien i den röda cirkeln så har kurvorna ovan parameterframställningen:

R    = 0.5
r    = R / 3         (3 = antalet utbuktningar)
c    = (R + r) / r   (= 3 + 1 = 4)
x(t) = c * r * cos(t) - q * r * cos(c * t)
y(t) = c * r * sin(t) - q * r * sin(c * t)   då 0 <= t <= 2 pi

Vi råkar ha en tjusig, 4.25 m lång, kantlist till bordet. Använd parameterframställningen ovan och bestäm q så att epicykloidens omkrets blir 4.25 m. Rita också upp formen på det resulterande bordet (glöm inte axis equal, annars får bordet fel form).

Ledning: Låt om(q) vara bordets omkrets (epicykloidens kurvlängd) som funktion av q. Då består problemet i att lösa ekvationen om(q) = 4.25. Använd fsolve för att lösa ekvationen där du använder integral för att beräkna kurvlängden.

Sist en liten övning som visar på för- och nackdelar med olika interpolationsförfaranden. Vi vill jämföra linjär interpolation, splineinterpolation och så kallad "shape-preserving" kubisk Hermiteinterpolation för två dataserier. Shape-preserving innebär att om data är monotont (växande/avtagande på ett delintervall) så är även interpolanten detta. En vanlig kubisk spline behöver inte alls vara shape-preserving.

De två dataserierna, (t(k), y1(k)), k = 1, ...,, 8, respektive (t(k), y2(k)), k = 1, ...,, 8, skapar du så här:

t = linspace(-2, 3, 8)'; y1 = exp(-t.^2); y2 = [-1 -1 -1 -1 1 1 1 1]';

I den första dataserien har du en underliggande funktion, exp(-t^2), att jämföra med när du studerar för- och nackdelar med de olika metoderna.

Använd Matlabrutinen interp1 för att testa de tre interpolationsmetoderna. Skapa ett plotfönster vardera för de två dataserierna. I varje plotfönster skall datapunkterna och de tre interpolaterna ritas ut. Datapunkterna skall ritas med ringar och interpolaterna skall ritas som kurvor på ett finare grid (använd 100 punkter). Använd slutligen legend-funktionen.

I den första dataserien vet du vilken funktion vi försöker approximera. Räkna ut maximala avvikelsen (på det finare gridet) för varje interpolant.

Slutligen: drag slutsatser och analysera dina resultat.