Integraler PDF

Preview

Summary

Integraler...

Description

Numerisk beräkning av integraler Trapetsregeln och Richardsonextrapolation Det finns en i det närmaste oöverskådlig mängd metoder för numerisk integralberäkning. Man kommer dock förvånansvärt långt med följande enkla betraktelse. Värdet av x2

integralen

f ( x)dx kan skattas med arean av en parallelltrapets (se fig.1 nedan) och vi x1

x2

får att

f ( x) dx  x1

f ( x1 )  f ( x 2 ) ( x2  x1 ) 2

fig.1

fig.2

Om denna approximation görs i en följd av intervall (se fig.2 ovan), där hela tiden x i 1  x i h får man xn

f ( x)dx  x0

f ( x 0 )  f ( x1 ) f ( x 1)  f ( x 2 ) f ( xn 1 )  f ( xn ) h h   h 2 2 2

f ( xn )   f ( x 0) h  f ( x1 )  f ( x2 )   f ( xn 1 )   2   2

Denna metod för approximativ beräkning av en bestämd integral har fått namnet Trapetsregeln. Vi betecknar Trapetsregeln med steget h med T (h) , alltså har vi f (x n )   f ( x0 )  f ( x 1 )  f ( x 2 )   f ( x n 1 )  T (h ) h   2   2 1

Exempel: Låt oss använda Trapetsregeln för att beräkna integralen

dx

1  x . 0

1 1 1 1 1 11      0.25 2.788095... 0.697024... Vi får T (0.25)   1  1.25 1.5 1.75 2 2  42

Integralens exakta värde kan beräknas till e log 2 0.693147... T (0.25) är som synes inte en särskilt bra skattning, men så har vi ju inte heller använt ett särskilt kort steg h .

Då steget h  0 gäller att T (h)  integralvärdet A . Trunkeringsfelet T (h )  A (om vi bortser från avrundningsfel) kan visas vara en potensserie i h av följande typ: T ( h)  A c 1h 2  c 2 h 4  c 3h 6  

Trunkeringsfelet vid användandet av Trapetsregeln är alltså av storleksordningen c1 h 2 (de övriga termerna är i stort sett försummbara i förhållande till denna term). Om vi byter ut h mot 2h i denna potensserieutveckling så får vi T ( 2h)  A c1 ( 2h) 2  c 2( 2 h) 4  c 3( 2h) 6   4c 1h 2 16c 2h 4  64c 3h 6  

Med hjälp av de två potensserieutvecklingarna kan vi nu eliminera bort termen med h 2 4(T (h)  A)  (T (2 h)  A)  4(c1 h2  c2 h4  c3 h6  )  (4c1 h2  16c 2 h4  64c3 h6  )  12c 2 h 4  60c 3 h 6  ...

vilket ger oss att T ( h) 

1 T (h )  T (2h)   A  4c 2 h 4  20c 3h 6   3

Så om vi beräknat T (h ) och T (2 h)så ger S ( h ) T ( h) 

1  T (h)  T (2h) en ännu 3

bättre approximation av integralens exakta värde A . Den teknik vi använt ovan för att eliminera feltermen c1 h 2 kallas lite mer allmänt för Richardsonextrapolation och kan användas i många andra sammanhang (än integralberäkning) för att förbättra numeriska metoder. Vi kan använda Richardsonextrapolation ytterligare gånger för att förbättra nogrannheten. I nästa steg elimineras termen med h 4 ur uttrycket för S (h )  A ovan (detta lämnas som övningsuppgift). Om man använder sådan upprepad Richardsonextrapolation på Trapetsregeln så kallas metoden för Rombergs metod. 1

Exempel: Låt oss använda uttrycket S (h) för att beräkna integralen

dx

1  x . 0

1 1 1 11    0.5 0.41666... 0.708333... och sedan tidigare Vi får T (0.5)   1  1.5 2 2  2 2 vet vi att T (0.25) 0.697024... så 1 S (0.25) T(0.25)   T (0.25)  T (0.5)  0.693253... 3

Detta är en förbättring av närmevärdet vi fick tidigare dvs. S (0.25) 0.693253... är en bättre approximation än T (0.25) 0.697024... av integralens exakta värde 0.693147...

Metoden med S (h ) innebär således att vi kan uppnå bättre nogrannhet på integralens värde utan att behöva beräkna funktionsvärdena i fler punkter än vad som krävs för att beräkna T (h ) . Låt oss undersöka S (h ) lite närmare. Använder vi att  f ( x0 ) f (xn )  T (h ) h   f (x1 )  f ( x2 )   f ( x n 1 )   så får vi 2   2 S( h) T( h)  

1  T( h)  T(2 h) 1  4T ( h)  T (2h)  3 3

h  f( x0)  4 f( x1)  2 f ( x2 )  4 f( x3 )    4 f( xn 3 )  2 f( xn 2 )  4 f ( xn 1 )  f ( xn) 3

Denna formel kallas Simpsons formel och är alltså den formel man får om man använder Richardsonextrapolation (en gång) på trapetsregeln. Man kan visa att samma formel dyker upp om man approximerar integranden med ett andragradpolynom, istället för förstagradspolynom som i Trapetsregeln.

Svårigheter vid numerisk integration Hur framgångsrika upprepade Richardsonextrapolationer kan bli vid numerisk integration sammanhänger med hur väl integranden lokalt kan approximeras med polynom. Det lönar sig att se efter om man med någon lämplig förbehandling av problemet kan få en integrand, som bättre lämpar sig för numerisk behandling. Låt oss titta på några exempel. Exempel (Serieutveckling): Antag att vi vill beräkna integralen

2

ex dx . Här kan inte  x 0

0.5

Trapetsregeln användas därför att integranden har en singularitet för x 0 . Om vi istället Maclaurinutvecklar integrandens täljare e

x2

  x 7 / 2 x 11 / 2 x 15 / 2 x 4 x6 x8  x  1 / 2 1  x 2     ...  x 1 / 2  x3 / 2     ... 2 6 24 2 6 24 x  

och integrerar denna serie term för term så får vi 0.5  x 2 e 2 2 2 2 0.517 /2 1.348143... Ju 0.59/ 2  0.513/ 2  dx 2 0.51/ 2  0.55/ 2   5 9 2 13 6 17 24    x 0 fler termer man tar med i serieutvecklingen, ju bättre approximation erhålls. Observera emellertid att serieutvecklingen av integranden endast ger bra närmevärden till integranden nära x 0 så om integrationsintervallet istället varit t.ex.  0,10  så skulle denna metod inte fungerat särskilt bra (såvida inte väldigt många termer i serieutvecklingen tagits med).

Exempel (Subraktion av singulariteten): Vi kan alternativt beräkna integralen i föregående exempel genom att från integranden subtrahera en funktion som har samma singularitet och som dessutom har en känd primitiv funktion. Till exempel har funktionen 1 / x samma singularitet som e  x / x och vi kan enkelt beräkna en primitiv till 1 / x . Lägger vi således till och drar ifrån denna funktion i integralen så får vi 2

0.5  x2

e



x

0

0.5

e

dx 

x2

0

 1 dx  x

0.5

 0

1 dx x

Den första av de två integralerna i högerledet ovan kan beräknas med t.ex. Trapetsmetoden, om vi definierar integrandens värde i x 0 som gränsvärdet 2

e x  1 0 . Den andra kan vi beräkna exakt. Vi har x 0 x

lim

0.5



1

0

1

Exempel (Substitution): Antag att vi vill beräkna integralen

x

ex

 0

dx  2 0.5 .

x

dx . Som i föregående

exempel kan vi inte heller här använda Trapetsregeln direkt eftersom integranden har en singularitet för x 0 . Om vi gör substitutionen x t 2 så försvinner emellertid 1

singulariteten ty

ex



x

0

1

et

dx   0

2

t

2

1

2tdt 2 et dt . Denna integral kan sedan beräknas 2

0

utan svårighet med Trapetsregeln (eller Simpsons formel om så önskas)

Exempel (Partiella integrationer): Vi kan ibland även använda partiell integration för att eliminera singulariteter. För integralen i föregående exempel har vi t.ex. 1

 0

ex x





1

1

1

0

0

dx  2 x e x 0  2 x e x dx  2e  2 x e xdx

Integralen i högerledet ovan har ingen singularitet så den kan beräknas med t.ex. Trapetsregeln. Om hög noggrannhet önskas så är det dock tillrådligt att göra ytterligare partiella integrationer, eftersom derivatan av den nya integranden har en singularitet i x 0 . Felet i Trapetsregeln beror ju på hur mycket funktionsvärdena varierar så om lutningen är stor blir approximationen sämre. Exempel (Kapning): Oändliga integrationsintervall förekommer ofta i praktiska problem. 

Antag att vi vill beräkna integralen e  x dx . Eftersom denna generaliserade integral är   konvergent så får vi en ganska god approximation om vi istället integrerar över något 2



begränsat intervall kring x 0 . Till exempel har vi att

4

e  x dx  e  x dx . Beräknar 2

2



4

vi den senare integral med Trapetsregeln för h 0.5 så får vi 1.772453... . Det exakta 

e x värdet på integralen 

2

dx

är

 1.772454...

 

Adaptiv integration Vårt mål är att med så få funktionsberäkningar som möjligt beräkna en given integral med föreskriven noggrannhet. Om integranden har stor andraderivata i en del av integrationsintervallet krävs en liten steglängd h för att få tillfredställande noggrannhet med trapetsregeln i detta delintervall, medan man i resten av intervallet kan få samma noggrannhet med mycket större h -värde. För att inte göra onödigt många funktionsberäkningar, måste man alltså anpassa steglängden h till integrandens b

c1

c2

a

a

c1

ck

uppförande. Vi gör således en uppdelning  ... 

b

   där integralerna över ck

1

ck

varje litet intervall behandlas som oberoende problem med olika val av steglängd. Låt oss titta på ett exempel. 10

Exempel: Antag att vi vill beräkna integralen

1

1 x

2

dx .

0

En första grov approximation av integralen får vi med Trapetsregeln för h 10 . 10

10  1 1   dx    5.049504... 2  1  0 2 1  10 2 

1

1  x

2

0

För att få en bättre approximation delar vi sedan in intervallet  0,10 i två lika stora delar, nämligen intervallen  0,5 och  5,10 , och beräknar integralerna över dessa delintervall med Trapetsregeln för h 5 . Vi får 5

5 1 1   dx     2.596153... 2 2  1 0 1 52 

1

1  x

2

0

10

1

1  x

2

5

5 1 1   dx    0.120906... 2 2  1 5 1  10 2 

Summerar vi dessa integraler så får vi att 10

1

1  x

2

dx 2.596153... 0.120906... 2.717060...

0

vilket naturligtvis är ett bättre närmevärde än den vi erhöll först. Men eftersom skillnaden av den första skattningen 5.049504... och den senare 2.717060 är relativt stor så behöver vi göra en finare indelning för att avgöra dess noggrannhet. Låt oss därför dela in var och en av intervallen  0,5 och  5,10 i två delintervall, nämligen  0,2.5 och  2.5,5 resp.  5,7 .5 och  7.5,10 , och beräkna integralerna över dessa delintervall med Trapetsregeln för h 2.5 . Vi får

2. 5

1

1  x

dx 

2

dx 

2

dx 

2.5  1 1    0.069910...  2 2 1 5 1  7.5 2 

2

dx 

2.5  1 1    0.034210...  2 2  1  7.5 1  10 2 

0

5

1

1 x

2. 5

7. 5

1

1  x

2.5  1 1     0.220490...  2 2  1  2.5 1 52 

5

10

2.5  1 1    1.422413...  2  1  0 2 1  2.5 2 

2

1

1  x

7. 5

Summerar vi de två första av dessa fyra integraler så får vi 5

1

1 x

2

dx 1.422413... 0.220490... 1.642904... (obs! tidigare fick vi

0

2.596153... ) och summerar vi de två övriga så får vi 10

1

1 x

2

dx  0.069910... 0.034210... 0.104121... (obs! tidigare fick vi

5

0.120906... )

Sammantaget får vi att 10

1 0 1 x 2 dx 1.642904... 0.104121... 1.747025.... (obs! tidigare fick vi 2.717060... )

Vi ser att närmevärdet till integralen över intervallet  0,5 avviker mycket mer från tidigare beräknat närmevärde av samma integral än den över intervallet  5,10 . Detta beror naturligtvis på att funktionen har mycket större andraderivata i intervallet  0,5 . Så för att få bättre noggrannhet nöjer vi oss med att dela in intervallen  0,2.5 och  2 .5,5 i ytterligare delintervall. Vi får 2.5

1.25 2.5 1 1 1 dx dx dx    2 2 0 1  x2   0 1 x 1.25 1  x

1.25  1 1 1 1   1.25  1.199011...      2 2  2 2  2  1 0 1  1.25  2  1  1.25 1 2.5  och 5 3.75 5 1 1 1    dx  2 dx 2 dx  x  2 x 3.75 1 x 2.5 1 2.5 1 



så

1.25  1 1 1 1   1.25        0.193232...   2 2 2 2  1  2.5 1  3.75  2  1  3.75 1 52 

5

1 1 x 2 dx 1.199011... 0.193232... 1.392244... (obs! tidigare fick vi 0  1.642904... ) 10

Tillsammans med tidigare beräknat värde av integralen

1

1  x

2

dx så får vi att

5

10

1

1 x

2

dx 1.392244 0.104121... 1.496365... (obs! tidigare fick vi

0

2.596153... )

Vi har fortfarande ganska dåligt närmevärde på integralen så ytterligare förfinad indelning av vissa delintervall är nödvändigt. Vi lämnar detta åt läsaren som förhoppningsvis har förstått idén bakom förfarandet ovan. Detta förfarande med en succesiv förfining av delintervall som ger dålig noggrannhet är exempel på en s.k. adaptiv metod. Följande Matlabkod implementerar en adaptiv integration med trapetsregeln. function I=trapets(f,a,b,tol,I0) m=(a+b)/2; h=b-a; fa=feval(f,a); fm=feval(f,m); fb=feval(f,b); I1=h*(fa+fm)/4; I2=h*(fm+fb)/4; I=I1+I2; if abs(I-I0)>tol; I1=trapets(f,a,m,0.5*tol,I1); I2=trapets(f,m,b,0.5*tol,I2); I=I1+I2; else hold on plot([a a m m],[0 f(a) f(m) 0]) plot([m b b],[f(m) f(b) 0]) x=a:h/100:b;plot(x,f(x),'r') hold off end

Om vi sparar ovanstående programrader i en m-fil med namnet trapets.m så kan t.ex. 10

integralen

1

1 x

2

dx beräknas med följande kommandon:

0

f=inline('1./(1+x.^2)');I=trapets(f,0,10,0.05,1) I = 1.4813 1

0. 9

0. 8

0. 7

0. 6

0. 5

0. 4

0. 3

0. 2

0. 1

0 0

1

2

3

4

5

6

7

8

9

10

Programmet är rekursivt (dvs. anropar sig självt) och integrationen på ett delintervall avbryts när skillnaden i två succesiva beräkningar av samma integral (I-I0) är tillräckligt liten (vilket bestäms av toleransen tol). Programmet trapets,m ger inte bara integralens approximativa värde, utan plottar även de trapetser vars areor beräknas. Observera att det krävs finare indelning på de delintervall där funktionsgrafen kröker mycket. I 1 6 x2  2 f x f x ( ) ''( )   och därmed . Genom att plotta ovanstående exempel är 1 x2 (1  x 2 ) 3 absolutbeloppet av andraderivatan får vi en uppfattning om var funktionsgrafen y  f ( x ) kröker mycket: x=0:0.01:10;df=(6*x.^2-2)./(1+x.^2).^3;plot(x,abs(df)) 2

1 . 8

1 . 6

1 . 4

1 . 2

1

0 . 8

0 . 6

0 . 4

0 . 2

0 0

1

2

3

4

5

6

7

8

9

1 0...