Kvadratur och interpolation


Matlabs quadl är en bra rutin som även kan hantera enklare ändpunkts-singulariteter. Om singulariteten blir för uttalad ger quadl ett felmeddelande och svaret blir mindre bra. Så här kan det se ut (jag har tagit bort varningen för nolldivision). Toleransen är 1e-6 absolut.

>> fel = quadl(@(x) 1 ./ x.^0.7, 0, 1) - 1 / 0.3, % 1 / 0.3 är integralens exakta värde
fel = -6.1687e-06                                 % analogt i följande exempel

>> fel = quadl(@(x) 1 ./ x.^0.8, 0, 1) - 1 / 0.2
Warning: Minimum step size reached; singularity possible.
fel =  -7.9354e-04

>> fel = quadl(@(x) 1 ./ x.^0.9, 0, 1) - 1 / 0.1
Warning: Minimum step size reached; singularity possible.
fel = -1.2600e-01

>> fel = quadl(@(x) 1 ./ x.^0.99, 0, 1) - 1 / 0.01
Warning: Minimum step size reached; singularity possible.
fel = -6.4572e+01

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 quadl 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 quadl för att approximera den bestämda integralen av cos(x) /  x^0.9.
b) Använd  variabelsubstitution följd av quadl 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 quadl 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.
 

Frivillig:
Kör (och studera) följande Matlabprogram. Titta på värdena och plotten. Förklara vad som händer.
Ledning: Hur lurar men en kvadraturmetod (vilken som helst)? Vikter och abscissor är givna (och fixa).
Skriv ett program för Simpson formel och beräkna en approximation av integralens värde.


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).

Epicyckloider
 
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 quadl 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"

Back