Opdracht Numerieke Methoden; Numerieke Integratie Juni 2011

Numerieke Methoden Numerieke Integratie

Den Helder, Juni 2011,

Cadet-sergeant Jeffrey L. A. Willems MST03 LVT OR

2

2

Samenvatting Dit paper is geschreven naar aanleiding van de opdracht ter numerieke berekening van integralen. Er worden verschillende methoden gebruikt om integralen numeriek o p te lossen. Een viertal methoden wordt in dit verslag uiteengezet. Zowel de trapeziumregel als een eigen ontworpen methode passeren de revue. Daarnaast wordt er ook gekeken naar standaardroutines van Matlab om deze integralen op te lossen. Deze zijn; QUAD en QUADL. De eigenmethode die ontworpen is, is net als de trapeziumregel een driepunts gesloten NewtonCotes formule, waarbij het verschil zit in de opdeling in het aantal deelintervallen. Dit heeft tot gevolg dat er een verschil zit in het aantal functie-evaluaties bij een gelijk blijvende stapgrootte h. Daarnaast is de ordebepaling afhankelijk van de continuïteit van een functie. Te zien is als een functie discontinu is dat de orde foutief bepaald wordt en de foutschatting afwijkt. De methoden die hiervoor gebruikt worden zijn dus ongeschikt voor discontinue functies. In tegenstelling tot continue functies, waarbij de orde en de foutschatting op juiste waarden bepaald worden.

3

Inhoud Samenvatting........................................................................................................... ........................................................................................................................................... ................................ 3 Inhoud .............................................................. ................................................................................................................................ ....................................................................................... ..................... 4 1.

Inleiding ............................................................... ................................................................................................................................. ............................................................................ .......... 6

2.

Numerieke Methoden ............................................................... ..................................................................................................................... ...................................................... 7 2.1

3.

Bepaling van gewichten uit algemene integratieformule ....................................................... 7

2.1.1

f(x)=1............................................................... ................................................................................................................................ ................................................................. 7

2.1.2

f(x)=x ................................................................ ................................................................................................................................ ................................................................ 8

2.1.3

f(x)=x .............................................................. ............................................................................................................................... ................................................................. 8

2.1.4

f(x)=x .............................................................. ............................................................................................................................... ................................................................. 8

2.1.5

Gewichtsbepaling middels lineaire algebra................................................. algebra..................................................................... .................... 9

2 3

2.2

Bepaling foutterm .............................................................. ................................................................................................................. ................................................... 10

2.3

Bepaling samengestelde integratieformule ........................................................ .......................................................................... .................. 11

2.4

Foutterm van samengestelde integratieformule. ......................................................... ................................................................. ........ 12

2.5

Orde van nauwkeurigheid en foutschatting................................................................. .......................................................................... ......... 13

2.6

Matlab implementatie................................................................. ........................................................................................................... .......................................... 14

2.7

Vergelijking met Trapeziumregel ........................................................... .......................................................................................... ............................... 14

2.8

Matlabfuncties ................................................................. ...................................................................................................................... ..................................................... 15

2.8.1

QUAD .............................................................. ............................................................................................................................. ............................................................... 15

2.8.2

QUADL ............................................................. ........................................................................................................................... .............................................................. 15

Experimenten .................................................................. ................................................................................................................................ .............................................................. 16 3.1

Functie 1; fx = 5 · x4 + 9 · x2 + 4 · x . ................................................................................. ................................................................................. 17

3.1.1

Numerieke integratie ............................................................. .................................................................................................... ....................................... 18

3.1.2

Numerieke ordebepaling van de methode .............................................................. ................................................................... ..... 19

3.1.2.1

QUAD ................................................................. ............................................................................................................................. ............................................................ 19

3.1.2.2

QUADL ................................................................ ........................................................................................................................... ........................................................... 20

3.2

Functie 2;

 

......................................................................................................... ................................................ 21 = sin7 4 .........................................................

3.2.1

PLOT g(x) ......................................................... ........................................................................................................................ ............................................................... 21

3.2.2

Numerieke integratie ........................................................... .................................................................................................... ......................................... 22

3.2.3

Standaard Matlabroutines .............................................................. ............................................................................................ .............................. 24

3.2.3.1

QUAD en QUADL op G(x) ................................................................. ............................................................................................... .............................. 24

3.2.3.2

QUAD en QUADL op H(x) ................................................................. ............................................................................................... .............................. 26

4

3.2.4

QUAD tegen QUADL op g(x) .......................................................................................... 28

4.

Conclusie ....................................................................................................................................... 29

5.

Literatuurlijst ................................................................................................................................. 30

Bijlage 1: Maple-uitwerking: opstellen vergelijkingen .......................................................................... 31 1.1

f(x)=1...................................................................................................................................... 31

1.2

f(x)=x ...................................................................................................................................... 31

1.3

f(x)=x^2 .................................................................................................................................. 32

1.4

f(x)=x^3 .................................................................................................................................. 32

Bijlage 2: Maple-uitwerking: bepaling coëfficiënten............................................................................. 33 Bijlage 3: Maple-uitwerking: bepaling hoogste macht .......................................................................... 34 Bijlage 4: Maple-uitwerking: hoogste afgeleide.................................................................................... 35 Bijlage 5: Maple-uitwerking: exacte functie-integralen ........................................................................ 36 Bijlage 6: Matlab-script: eigenmethode ................................................................................................ 37 Bijlage 7: Matlab-script: Trapeziummethode........................................................................................ 38 Bijlage 8: Matlab-script: Integraal 1 ...................................................................................................... 39 Bijlage 9: Matlab-script: Integraal 2 ...................................................................................................... 42 Bijlage 10: Matlab-script: Integraal 3 .................................................................................................... 45

5

1. Inleiding In dit paper worden verschillende numerieke integratietechnieken uiteengezet. Er worden praktische voorbeelden uitgewerkt en verschillende technieken toegepast op probleemstellingen. Deze technieken hebben de overeenkomst de oppervlakte onder een bepaalde grafiek te berekenen. Dit aan de hand van gestandaardiseerde vormen, afhankelijk van de functiewaarden en lengte van het te integreren oppervlak. Ook wordt er gebruik gemaakt van computeranalyses met behulp van wiskundige software. In dit verslag wordt het softwareprogramma MATLAB gebruikt. Hierin zitten standaard functies die numeriek integralen op kunnen lossen. Daarnaast wordt in dit verslag methoden tegen elkaar uitgezet en vergeleken. Zo kan de implementatie van deze functies met elkaar vergeleken worden en hieruit een o ptimalisatie kiezen. In hoofdstuk 2 staat de theorie uiteengezet, inclusief de uitleg hoe de standaard matlab routines werken. In hoofdstuk 3 wordt de theorie toegepast een drietal functies. In het laatste hoofdstuk worden conclusies getrokken uit de verkregen resultaten.

6

2. Numerieke Methoden Numerieke Integratie In dit hoofdstuk wordt een algemene integratieformule uitgewerkt (2.1). van deze integratieformule worden de gewichten A, B, C en D bepaald. Naast deze gewichten wordt ook de foutterm uitgewerkt. Binnen deze foutterm worden de waarden van K, p en s bepaald. De foutterm is gelijk aan:

 ℎ+1    ·

·

( )

( )

De eis die geldt voor parameters binnen de integratieformule is dat de functie geldt voor polynomen met een zo hoog mogelijke graad. Verder wordt in dit hoofdstuk de samengestelde integratieformule inclusief de foutterm op basis van N deelintervallen van lengte 4h bepaald. Ook wordt de orde van de nauwkeurigheid en een formule van de foutschatting bepaald. Achtereenvolgens worden de gevonden integratieformules in Matlab geïmplementeerd door een correctie op een bestaande matlab-file.

2.1

Bepaling van gewichten uit algemene integratieformule

Om de gewichten uit formule 1 te bepalen moet de integraal gelijk gesteld worden aan de algemene formule. Door te variëren in de functie f(x) binnen de integraal kan een stelsel vergelijkingen opgesteld worden, waaruit door middel van lineaire algebra de waarden van A, B, C en D opgelost kunnen worden. Hierbij moet de foutterm gelijk zijn aan nul.

∫−22       ℎ   ℎ     ℎ  ℎ+1    ( )

=

· ( 2 )+

· (

)+

· (0) +

· (2 ) +

·

·

( )

( )

(2.1)

2.1.1 f(x)=1 De eerste vergelijking die opgesteld wordt is door de algemene formule gelijk te stellen aan de integraal van -2h tot 2h over de integraal van f1(x)=1. Daarna wordt deze integraal gelijk gesteld aan de algemene formule. Voor de uitwerkingen in Maple wordt verwezen naar bijlage 1.

2 ∫−2    ℎ    ℎ   ℎ     ℎ ℎ     ℎ 1( ) = 1

1( )

=4

· ( 2 )+ ·1+

·1+

· (

·1+

7

)+

· (0) +

·1=4

· (2 ) = 4

(2.2) (2.3) (2.4) (2.5)

2.1.2

f(x)=x

In plaats van f(x)=1 wordt de vergelijking van f(x)=x geïmplementeerd in de algemene vergelijking

2  ∫−2       ℎ   ℎ     ℎ   ℎ  ℎ  ℎ 2( ) =

2( )

2.1.3

=0

· ( 2 )+

· (

2 ·

+2

·

)+

·

· (0) +

· (2 ) = 0

=0

(2.7) (2.8) (2.9)

f(x)=x2

2 2 16 ∫−2    3 ℎ3    ℎ   ℎ     ℎ 163 ℎ3  ℎ2  ℎ2  ℎ2 163 ℎ3 3( ) =

3( )

=

· ( 2 )+

4 ·

2.1.4

(2.6)

+

·

· (

+4

)+

·

· (0) +

· (2 ) =

=

(2.10) (2.11) (2.12) (2.13)

f(x)=x3

2 3 ∫−2      3 ℎ  3 ℎ 3    ℎ   ℎ  ℎ  ℎ 4( ) =

3( )

=0

· ( 2 )+

8 ·

·

· (

)+

+8 ·

8

· (0) +

=0

· (2 ) = 0

(2.14) (2.15) (2.16) (2.17)

2.1.5

Gewichtsbepaling middels lineaire algebra

Via het verkregen stelsel vergelijkingen kunnen door middel van lineaire algebra de waarden van A, B, C en D berekend worden. Door het stelsel in matrixvorm te zetten en deze matrix terug te brengen naar een reduced echelon vorm, kunnen de waarden uitgelezen worden. Voor de mapleberekeningen wordt verwezen naar bijlage 2.

ℎ ℎℎ2ℎ3 ℎℎ2ℎ3 ℎℎℎ23 ⎣ ℎ3⎦ ℎ ℎ ℎ⎦ ⎣ 1 2 4 8

Reduced row echalon form

1 1 1 1

1 0 0 0

1

0

0

0

0

1

0

0

0

0

1

0

0

0

0

1

Hiermee kunnen de waarden uitgelezen worden:

 23 ℎ 8  3ℎ  23 ℎ = · =0

= ·

= ·

9

4 0 16

1 2 4 8

2 3

3

0

·

0

8 3 2 3

· ·

2.2

Bepaling foutterm

Om de foutterm te bepalen wordt gekeken naar de orde n van de functie f(x);

   ( )=

waarbij de integraal van de functie f(x) niet gelijk is aan de algemene vergelijking.

  4 2 4 64 ∫−2    5 ℎ5 23 ℎ   ℎ 83 ℎ  23 ℎ  ℎ 643 ℎ5 ∫−22       ℎ  64 ℎ 64     ℎ  ℎ+1    5 ℎ5 3 ℎ5  ℎ+1    Na volgende berekeningen is te zien dat bij

( )=

de integraal niet gelijk is aan algemene

formule. Voor maple-berekeningen wordt verwezen naar bijlage 3.

( )=

( )

=

(2.18) (2.19)

Substitutie van coëfficiënten in formule 2.1 geeft:

·

· ( 2 )+ ·

· (0) + ·

· (2 ) =

(2.20)

Teruggebracht naar vergelijking 2.1 geeft:

( )

=

· ( 2 )+

· (

)+

=

· (0) + +

·

· (2 ) + ·

( )

·

·

( )

( )

( )

(2.1) (2.21)

Om geen afhankelijke functiewaarde van h in het rechter lid van de vergelijking te krijgen, moet gekeken worden naar de hoogste afgeleidde van de functie f(x). Om geen afhankelijke waarde van h te behouden is de hoogste afgeleide en waarde van ‘s’ gelijk aan 4. Voor de berekening, zie bijlage 4. (4)

De waarde van f (c) is gelijk aan 24. Om de waarde van ‘K’ te bepalen moet gekeken worden naar formule 2.21. Hiervoor moet eerst de waarde van ‘P’ bepaald worden. Doordat ‘K’ de correctiefactor is, moet de hoogste macht van ‘h’ gelijk zijn aan 5, waarop de waarde van ‘p’ op 4 uitkomt. Hiermee komt vergelijking 2.21 op het volgende uit:

645 ℎ5 645 ℎ5

643 ℎ5  ℎ4+1 643 ℎ5  ℎ5      −     24  1645 83 ℎ  23 ℎ  ℎ  1645 ℎ5  4  =

+

·

=

+ 24

· 24

(2.22)

·

De waarde van K wordt volgens formule 2.23 berekend:

=

·

=

(2.23)

Hiermee komt de algemene formule uit op (formule 2.1+2.20+2.23):

∫−22    23 ℎ   ℎ ( )

= ·

· ( 2 )+ ·

· (0) + ·

· (2 )

10

·

·

( )

( )

(2.24)

2.3

Bepaling samengestelde integratieformule

De samengestelde integratieformule inclusief de foutterm wordt op basis van N deelintervallen bepaald, waarbij de deelintervallen een lengte van 4h hebben. Er wordt geïntegreerd over het interval [a,b] over N deelintervallen. Voor de benadering van lengte 4h. Kortom:

∫     ℎ    ℎ    ( )

wordt het interval over [a,b] dus over N deelintervallen van

·4 =

,

·

=

4

·

Hierbij loopt K van 1 tot N over N deelintervallen.

De samengestelde integratieformule wordt dan:

∫−22    ℎ23   83   ℎ 23   ℎ  1645 ℎ5  4  ( )

=

· ( )+ · ( +2 )+ · ( +4 )

·

·

·

( )

(2.25)

( )

Met de deelintervallen geïntegreerd geeft:

∫−22    ∑=1 ℎ23  ℎ    83  ℎ    23   ℎ  1645 ℎ5  4  ( )

=

+4 (

·

1) + ·

+ (4

2)

+ · ( +4

)

·

·

·

( )

(2.26)

( )

Te zien is dat het interval over 4h verdeeld wordt, terwijl er op (a+h) en (a+3h) de gewichten nul gelden. Hierdoor kan een andere samengestelde integratieformule opgesteld worden, waardoor het aantal functie-evaluaties afneemt.

 2 −2   =1 ℎ�   ( )

=

2 3

· ( )+

 ℎ      ℎ    ℎ5  4 

8 3

·

+ (4

2)

11

+

4 3

· ( +4

)+

2 3

· ( )

·

16 45

·

·

( )

( )

Deze formule kan vereenvoudigd worden door het aantal functie-evaluaties terug te brengen naar 2N+1 bij N intervallen:

 2 −2   =1 ℎ�   ( )

2.4

2

=

3

· ( )+

 ℎ      ℎ    ℎ5  4 

8 3

+ 2 (2

·

1)

+

4 3

· ( +4

)+

2 3

· ( )

·

16 45

·

·

( )

( )

Foutterm van samengestelde integratieformule.

De foutterm van de samengestelde formule komt uit o p:

∑=1  ℎ    ℎ5   4   ℎ    ℎ    ℎ4 −4  4   454 ℎ4    4  16 45

·

5

·

(4)

( )

=

16 45

·

·

·

( )

(2.27)

( )

Waarbij geldt dat;

·4 =

,

·

=

4

Gesubstitueerd geldt: 16 45

·

·

·

( )

( )

(2.28)

Wat te vereenvoudigen valt naar:

·

·(

)·

( )

( )

4

Hieruit valt op te maken dat bij een halvering van h, de fout 2 , 16 maal kleiner wordt.

12

(2.29)

2.5

Orde van nauwkeurigheid en foutschatting

Om de orde van de functie te bepalen wordt stap voor stap h verkleind. Voor een methode M van de p

orde h geldt:

  ℎ ≈  ℎ −  −  ≈  ( )

(2.30)

·

Hieruit volgt dat ( ) ( )

(2.31)

2

Echter kan de exacte waarde van W nier bepaald worden, daardoor wordt de formule omgeschreven naar waarden die wel te bepalen zijn, namelijk waarden M(2h), M(h/2) en M(h):

 ℎℎ  ℎℎ ( )

(2 )

=2

( )

2



Uit deze formule kan de waarde van p berekend worden

2

log

 ℎℎ   ℎℎ ( )

(2 )

=

( )

2

De foutschatting van de formule kan door de Richardson-extrapolatie [1, p. 181] bepaald worden. Voor deze methode moet echter de orde van de integratiemethode bekend zijn.

h W M ( ) 2 −

≈

1

h h 1 ( M ( ) − M ( h)) ≈ ( M ( ) − M ( h)) 2 −1 2 15 2 p

13

(2.32)

2.6

Matlab implementatie

Voor het matlabscript, betreffende de eigenmethode, wordt verwezen naar bijlage 6. Gegeven is dat het interval van N over 4h verdeeld wordt, dus:

ℎ   =

hierbij loopt k van 0 tot 4*N

2.7

(

)

4

Vergelijking met Trapeziumregel

De samengestelde trapeziumregel behoort tot de klasse van gesloten Newton-Cotes formules ([1] p. 191), dit vanwege de gelijke afstanden tussen de waarden van h onderling en omdat de eindpunten 2

behoren tot de punten die berekend worden voor de integratie. Deze methode heeft een fout O(h ) waardoor geldt dat de trapeziumregel exact is voor lineaire functies, dit wil zeggen dat de regel van de tweede orde is. Daarnaast heeft de methode voor N intervallen N+1 functie-evaluaties nodig. Voor de ontworpen integratie formule, ofwel de eigenmethode, zijn er 2N+1 functie-evaluaties nodig, maar deze is wel exact voor polynomen tot de graad vier. Een overeenkomst tussen beide methoden is dat het gesloten driepunts Newton-Cotes formules zijn. Door het verschil in de orde van beide functies zal de foutschatting van de zelf ontworpen methode kleiner zijn dan de foutschatting van de trapeziummethode; Foutschatting trapeziummethode:

13 12 ℎ ℎ  115 12 ℎ ℎ  ( )

( )

Foutschatting ontworpen methode:

14

2.8

Matlabfuncties

In matlab zitten voorgeprogrammeerde functies om numeriek integralen op te lossen. Twee van deze functies zijn QUAD en QUADL.

2.8.1 QUAD QUAD, wat staat voor quadrature (kwadratuur) maakt adaptief gebruik van de Simpsonregel. Deze aanpassing houdt in dat de functie QUAD de grootte van de intervallen af laat hangen van de richtingscoëfficiënt van de te integreren functie. Dit kan ook handmatig ingevuld worden, door bij de functie meerdere waarden in te voeren waartussen geintergreerd moet worden. Ook is het mogelijk -6

om de tolerantie op te geven. Standaard staat deze ingesteld op 1e . Daarnaast is de mogelijkheid om het aantal functie-evaluaties aanwezig door de TRACE-optie t e gebruiken. [Q,FCNT] = QUAD(FUN,A,B,TOL,TRACE) Het nadeel van de QUAD-functie is dat bij een wisselende functie met hoge afgeleiden de nauwkeurigheid van de integraal afneemt. Een oplossing is om in plaats van QUAD, QUADL te gebruiken.

2.8.2 QUADL QUADL berekent numeriek te integraal van de ingevoerde functie middels adaptieve Lobatto kwadratuur. De functie zoekt binnen de ingevoerde grenzen van de functie grote veranderingen en neemt deze waarden als tussenwaarden waarbinnen geïntegreerd wordt. [Q,FCNT] = QUADL(FUN,A,B,TOL,TRACE)

15

3. Experimenten Toepassing van theorie op praktijkvoorbeelden In dit hoofdstuk worden de methoden die in de voorgaande hoofdstukken zijn besproken geïmplementeerd. Deze methoden worden middels MATLAB gebruikt om een drietal integralen op te lossen. Bij het oplossen worden de stapgroottes van h gevarieerd. Tevens zal, zoals besproken in hoofdstuk 2, de nauwkeurigheid, de orde en foutschatting berekend worden. Tevens worden de functies QUAD en QUADL gebruikt om de integralen op te lossen.

4 2 1       5 3 2 7 4   √  ∫01 √ 7 4   ∫12 √ 7 4  In paragraaf 3.1 zal de functie ( ) = 5 · integraal F:

=

+3·

+2·

+9·

+4·

behandeld worden, met bijbehorende

.

Omdat de ontworpen methode nauwkeurig is tot en met de derde graad zal er een fout ontstaan. Deze fout kan eenvoudig bepaald worden om dat de integraal van de functie exact te berekenen is. In paragraaf 3.2 zal de functie

=

sin

( ) = sin

behandeld worden, met bijbehorende integraal G:

en in paragraaf 3.3 de integraal H:

16

=

sin

3.1

4

2

Functie 1; f (x) = 5 · x + 9 · x + 4 · x .

De volgende functie wordt als eerste experiment genomen;

4 2 ∫02    ∫02 4 2

f (x) = 5 · x + 9 · x + 4 · x

(3.1)

de integraal die hierbij opgelost moet worden is gelijk aan:



=

( )

=

5 3 2 20

5 · x + 9 · x + 4 · x dx = [

+3·

+2·

De integraal is grafisch weergegeven in figuur 1, zie onderstaand.

Figuur 1: Plot integrand F(x)

17

] = 64

(3.2)

3.1.1

Numerieke integratie

De integraal F wordt tevens berekend met behulp van MAPLE in Bijlage 5. Deze waarde wordt gebruikt als de exacte waarde om de ware fout te berekenen. De gegevens verkregen met MATLAB zijn in de volgende tabellen weergegeven met behulp van het script wat in bijlage 8 is bijgevoegd. Tabel 1: Gegevens trapeziumregel van functie f(x)

h

Nf

Benadering

Geschatte fout

Ware fout

Orde

1/8

17

64,255126953125

-0,254801432292

0,255126953125

1,998273117488

1/16

33

64,063796997070

-0,063776652018

0,063796997070

1,999568602155

1/32

65

64,015950202942

-0,015948931376

0,015950202942

1,999892170692

1/64

129

64,003987610340

-0,003987530867

0,003987610340

1,999973043932

1/128

257

64,000996906310

-0,000996901343

0,000996906310

1,999993261062

Geschatte fout

Ware fout

Orde

Tabel 2: Gegevens eigenmethode van functie f(x)

h

Nf

Benadering

1/8

9

64,000325520833

1/16

17

64,000020345052

-0,000020345052

0,000020345052

4,000000000000

1/32

33

64,000001271566

-0,000001271566

0,000001271566

4,000000000000

1/64

65

64,000000079473

-0,000000079473

0,000000079473

4,000000000000

1/128

129

64,000000004967

-0,000000004967

0,000000004967

4,000000000000

0,000325520833

Ten eerste is op te merken dat voor eenzelfde stapgrootte (beginnend met h/8) er minder functieevaluaties nodig zijn als men gebruik maakt van de eigenmethode. Dit is vanwege dat de eigenmethode per N deelinterval het deelinterval opdeelt in 4h, terwijl de trapeziummethode dit niet doet. Hierdoor moet N dan ook verschillend worden gekozen per methode (zie bijlage 8 voor de verschillen in N deelintervallen), namelijk met een verdeling van 1:4. Daarnaast is te zien dat de eigenmethode een meer nauwkeurige benadering geeft ten opzichte van de trapeziumregel en dat de geschatte fout niet afwijkt van de ware fout, hier is uit op te maken dat de methode die gebruikt wordt (zie formule 2.32) voor de foutschatting een nauwkeurige benadering geeft van de ware fout.

18

3.1.2 Numerieke ordebepaling van de methode Zie kolommen 6 uit voorgaande tabellen. Te zien is dat de orde van de trapeziummethode gelijk is aan 2 en de eigenmethode gelijk is aan 4, wat overeenkomt met de theorie behandeld in paragraaf 2.7.

3.1.2.1

QUAD

Tabel 3: QUAD op F(x)

Tolerantie

Benadering

Fout

Nf

1,00E-01

64,000000000000

0,000000000000

13

1,00E-02

64,000000000000

0,000000000000

13

1,00E-03

64,000000000000

0,000000000000

17

1,00E-04

64,000000000000

0,000000000000

25

1,00E-05

64,000000000000

0,000000000000

33

1,00E-06

64,000000000000

0,000000000000

65

1,00E-07

64,000000000000

0,000000000000

97

1,00E-08

64,000000000000

0,000000000000

129

1,00E-09

64,000000000000

0,000000000000

257

1,00E-10

64,000000000000

0,000000000000

385

1,00E-11

64,000000000000

0,000000000000

513

1,00E-12

64,000000000000

0,000000000000

1025

1,00E-13

64,000000000000

0,000000000000

1537

1,00E-14

64,000000000000

0,000000000000

2049

1,00E-15

64,000000000000

0,000000000000

4097

Uit voorgaande tabel is op te maken dat de functie QUAD de integraal exact bepaald met het aantal functie-evaluaties toenemend met een aflopende tolerantie.

19

3.1.2.2

QUADL

Tabel 4:QUADL op F(x)

Tolerantie

Benadering

Fout

Nf

1,00E-01

64,000000000000

0,000000000000

18

1,00E-02

64,000000000000

0,000000000000

18

1,00E-03

64,000000000000

0,000000000000

18

1,00E-04

64,000000000000

0,000000000000

18

1,00E-05

64,000000000000

0,000000000000

18

1,00E-06

64,000000000000

0,000000000000

18

1,00E-07

64,000000000000

0,000000000000

18

1,00E-08

64,000000000000

0,000000000000

18

1,00E-09

64,000000000000

0,000000000000

18

1,00E-10

64,000000000000

0,000000000000

18

1,00E-11

64,000000000000

0,000000000000

18

1,00E-12

64,000000000000

0,000000000000

18

1,00E-13

64,000000000000

0,000000000000

18

1,00E-14

64,000000000000

0,000000000000

18

1,00E-15

64,000000000000

0,000000000000

108

Net als de functie QUAD bepaalt QUADL de integraal van de polynoom tussen [0 2] de waarde exact. Wel is op te merken dat QUADL minder functie-evaluaties nodig heeft (bij lagere tolerantie).

20

3.2

  √ 7 4

Functie 2; ( ) = sin

De volgende functie wordt als tweede experiment genomen;

  7 4 1 7 4  0     ( ) = sin

De eerste integraal die hierbij opgelost moet worden is gelijk aan:

=

sin

De tweede integraal die hierbij opgelost moet worden is gelijk aan:

 3.2.1

=

2 7 4 1     sin

PLOT g(x)

Figuur 2: plot Functie g(x)

21

3.2.2

Numerieke integratie

In de volgende tabellen zijn de gegevens betreffende de verschillende methoden uiteengezet. Net als bij de functie f(x)is er een verschil in het aantal deelintervallen opdat eenzelfde grootte van h te behouden. Hierdoor verschilt ook het aantal functie-evaluaties. Wat wel op te merken is, is dat de orde van zowel de trapeziummethode als de eigenmethode toegepast op de functie G(x) afwijkt van de theoretische waarde. Waar dit aan zou kunnen liggen is de functie waarover geïntegreerd wordt discontinue is op [0,0]. Als dan gekeken wordt naar de tabellen 7 en 8 is te zien dat de ordes wel weer overeenkomen met de theoretische waarden van twee en vier. Naast het verschil in de orde zit er binnen de methoden ook een verschil in de ware fout en de geschatte fout. Ook dit is te wijten aan de discontinuïteit van de functie. Dit is eveneens, net als de orde, niet te zien als geïntegreerd wordt over de interval [1,2]. Hiermee valt te concluderen dat de afwijking ligt aan de discontinuïteit van de functie op het interval [-0,1 0,1] rond het punt (0,0)

3.2.2.1

Methoden op G(x)

Tabel 5: Trapeziumregel op G(x)

h

Nf

Benadering

Geschatte fout

Ware fout

Orde

1/8

9

0,570497516672

1/16

17

0,574830073252

0,001444185527

0,002251310817

1,544835437993

1/32

33

0,576314991136

0,000494972628

0,000766392933

1,552265665193

1/64

65

0,576821309770

0,000168772878

0,000260074299

1,557485209728

1/128

129

0,576993327870

0,000057339367

0,00008805619 9

1,561215281272

1/256

257

0,577051618871

0,000019430334

0,00002976519 8

1,563913864170

Geschatte fout

Ware fout

Orde

0,006583867397

Tabel 6: Eigenmethode op G(x)

h

Nf

Benadering

1/8

5

0,56992905711679 6

1/16

9

0,57467979448129 6

0,0003167158243000 25

0,00240158958770376

1,57507965883746

1/32

17

0,57627425877842 1

0,0001062976198082 93

0,00080712529057935 8

1,57356116681605

1/64

33

0,57680996376357 6

3,57136656770211e-05

0,00027142030542404 2

1,57248927303280

1/128

65

0,57699008264813 6

1,20079256373066e-05

9,13014208644425e-05

1,57192670781855

1/256

129

0,57705066723612 2

4,03897253245707e-06

3,07168328775864e-05

1,57165743215592

0,00715232695220414

22

3.2.2.2

Methoden op H(x)

Tabel 7: Trapeziumregel op H(x)

h

Nf

Benadering

Geschatte

fout

Ware

fout

(1.0e-

Orde

1/8

9

0.940927539677032

1/16

17

0.941193825667812

0.887619969266821

0.088789602188055

1,999579508824

1/32

33

0.941260416571415

0.221969678676492

0.022198698585107

1,999894553547

1/64

65

0.941277065514141

0.055496475753388

0.005549755859091

1,999973617974

1/128

129

0.941281227825936

0.013874372651509

0.001387444063639

1,999993403005

1/256

257

0.941282268408643

0.003468609023708

0.000346861356526

1,999998351160

0.355075592968102

Tabel 8: Eigenmethode op H(x)

h

Nf

Benadering

1/8

5

0,941275933318869

1/16

9

0,941282178884440

1/32

17

1/64

Geschatte

fout

Ware fout (1.0e-007)

Orde

6,68195113129322e-06

3,93343465076852

4,16371038092223 e-07

4,36385559909880e-07

3,98172033853617

0,941282587664739

2,72520198999852 e-08

2,76052614101019e-08

3,99530165863411

33

0,941282613539282

1,72496958198802 e-09

1,73071768028166e-09

3,99882084260850

1/128

65

0,941282615161716

1,08162271731752 e-10

1,08283604305370e-10

3,99964650966654

1/256

129

0,941282615263201

6,76566950611838e-12

6,79856171359461e-12

3,93343465076852

23

3.2.3 Standaard Matlabroutines In deze paragraaf worden de MATLAB routines weergegeven en besproken. Net als in paragraaf 3.1 is er gekozen voor aflopende tolerantie van 1e-1 tot 1e-15. Hierbij wordt gekeken naar de benadering van de integraal, de fout (exacte waarde berekend via MAPLE minus de benadering) en het aantal functie-evaluaties nodig om tot een waarde te komen die binnen de tolerantie valt. Verder is een grafiek weergegeven die de verschillen in de fout ten opzichte van het aantal functie-evaluaties weergeeft.

3.2.3.1

QUAD en QUADL op G(x)

Tabel 9: QUAD op G(x)

Tolerantie

Benadering

Fout

Nf

1,00E-01

0,576281624320449

0,000799759749

13

1,00E-02

0,576281624320449

0,000799759749

13

1,00E-03

0,576281624320449

0,000799759749

13

1,00E-04

0,576811312615695

0,000270071453

17

1,00E-05

0,577049343770596

0,000032040298

25

1,00E-06

0,577077840888460

0,000003543181

37

1,00E-07

0,577080969250011

0,000000414819

53

1,00E-08

0,577081337126668

0,000000046942

85

1,00E-09

0,577081378744645

0,000000005324

137

1,00E-10

0,577081383480327

0,000000000589

209

1,00E-11

0,577081384000814

0,000000000068

333

1,00E-12

0,577081384064450

0,000000000005

541

1,00E-13

0,577081384066849

0,000000000002

817

1,00E-14

0,577081384067121

0,000000000002

1305

1,00E-15

0,577081384067152

0,000000000002

2129

24

Tabel 10: QUADL op G(x)

Tolerantie

Benadering

Fout

Nf

1,00E-01

0,575986837343447

0,001094546726

18

1,00E-02

0,575986837343447

0,001094546726

18

1,00E-03

0,575986837343447

0,001094546726

18

1,00E-04

0,577055759722241

0,000025624347

48

1,00E-05

0,577080783264405

0,000000600805

78

1,00E-06

0,577080783264405

0,000000600805

78

1,00E-07

0,577081369987185

0,000000014082

138

1,00E-08

0,577081383726471

0,000000000343

168

1,00E-09

0,577081383737205

0,000000000332

228

1,00E-10

0,577081384059210

0,000000000010

318

1,00E-11

0,577081384059424

0,000000000010

408

1,00E-12

0,577081384066974

0,000000000002

588

1,00E-13

0,577081384067151

0,000000000002

888

1,00E-14

0,577081384067151

0,000000000002

1128

1,00E-15

0,577081384067155

0,000000000002

1548

Uit de gegevens van de tabellen 9 en 10 kunnen de ware fout ten opzichte van het aantal functieevaluaties tegen elkaar uitgezet worden, de grafiek die hieruit volgt is in de volgende figuur weergegeven. Hierbij stelt de nauwkeurigheid de ware fout voor.

Nauwkeurigheid tegen f-evaluaties QUAD en QUADL

-2

10

quad quadl

-3

10

-4

10

-5

10 ) t u o f

e r a w (

-6

10

d -7 i e 10 h g i r u e -8 k w 10 u a n -9

10

-10

10

-11

10

-12

10

1

10

2

3

10

4

10

10

f-evaluaties

Figuur 3: f-evaluaties tegen ware fout op G(x)

25

3.2.3.2

QUAD en QUADL op H(x)

Tabel 11: QUAD op H(x)

Tolerantie

Benadering

Fout (1.0e-009)

Nf

1,00E-01

0,941282611474996

3,795004421825

13

1,00E-02

0,941282611474996

3,795004421825

13

1,00E-03

0,941282611474996

3,795004421825

13

1,00E-04

0,941282611474996

3,795004421825

13

1,00E-05

0,941282611474996

3,795004421825

13

1,00E-06

0,941282611474996

3,795004421825

13

1,00E-07

0,941282614796352

0,473648342769

17

1,00E-08

0,941282615179442

0,090558116561

21

1,00E-09

0,941282615262245

0,007754796805

33

1,00E-10

0,941282615269574

0,004261035969

57

1,00E-11

0,941282615269946

0,000053956839

81

1,00E-12

0,941282615269967

0,000033306691

129

1,00E-13

0,941282615269969

0,000031308289

225

1,00E-14

0,941282615269969

0,000031308289

321

1,00E-15

0,941282615269969

0,000031308289

513

Tolerantie

Benadering

Fout (1.0e-011)

Nf

1,00E-01

0,941282615194927

7,507305888055

18

1,00E-02

0,941282615194927

7,507305888055

18

1,00E-03

0,941282615194927

7,507305888055

18

1,00E-04

0,941282615194927

7,507305888055

18

1,00E-05

0,941282615194927

7,507305888055

18

1,00E-06

0,941282615194927

7,507305888055

18

1,00E-07

0,941282615194927

7,507305888055

18

1,00E-08

0,941282615194927

7,507305888055

18

1,00E-09

0,941282615194927

7,507305888055

18

1,00E-10

0,941282615194927

7,507305888055

18

1,00E-11

0,941282615269969

0,003141931160

48

1,00E-12

0,941282615269969

0,003141931160

48

1,00E-13

0,941282615269969

0,003141931160

48

1,00E-14

0,941282615269969

0,003141931160

48

1,00E-15

0,941282615269969

0,003130828929

138

Tabel 12: QUADL op H(x)

26

Net als op de functie G(x) kan op het interval van H(x) een plot gemaakt worden van het aantal functie-evaluaties tegen de ware fout. De gegevens zijn in het volgende figuur weergegeven:

Nauwkeurigheid tegen f-evaluaties QUAD en QUADL

-8

10

quad quadl -9

10

-10

10

) t u o f e r a w ( d -11 i e 10 h g i r u e k w u a n -12

10

-13

10

-14

10

1

10

2

3

10 f-evaluaties

10

Figuur 4: f-evaluaties tegen ware fout op H(x)

27

3.2.4

QUAD tegen QUADL op g(x)

In figuur drie is te zien dat de functie QUAD een lagere ware fout heeft bij een lager aantal functieevaluaties. Naarmate te evaluaties toenemen, neemt de ware fout van zowel QUAD als QUADL af, maar de sterkte van de toename van QUADL is hoger. Hiermee valt te concluderen dat voor een nauwkeurigere waarde van de integraal er m eer f-evaluaties nodig zijn. Ook is uit de figuren op te maken dat de functie QUAD beter werkt voor discontinue functies, dit aan de lagere foute waarde op kleiner aantal evaluaties. Voor continue functies moet er gebruik gemaakt worden van de functie QUADL.

28

4. Conclusie Dit paper is geschreven om inzicht te krijgen in het numeriek oplossen van integralen. Deze functies kunnen zowel continuïteit als discontinuïteit bevatten. In hoofdstuk twee staat de theorie beschreven omtrent de orde van een numerieke integratiemethoden en de foutschatting. De eigenmethode is net als de trapeziumregel een driepunts gesloten Newton-Cotes formule, waarbij het verschil zit in de opdeling in het aantal deelintervallen. Dit heeft tot gevolg dat er een verschil zit in het aantal functie-evaluaties en deelintervallen bij een gelijk blijvende stapgrootte h. Daarnaast is de ordebepaling afhankelijk van de continuïteit van een functie. Te zien is als een functie discontinu is dat de orde foutief bepaald wordt en de foutschatting afwijkt. De methoden die hiervoor gebruikt worden zijn dus ongeschikt voor discontinue functies. In tegenstelling tot continue functies, waarbij de orde en de foutschatting op juiste waarden bepaald worden.

29

5. Literatuurlijst th [1] Burden, R.L., Faires, J.D. , Numerical analysis 8 edition, Thomson Brooks/Cole, Belmont, 2005.

30

Bijlage 1: Maple-uitwerking: opstellen vergelijkingen

1.1

f(x)=1

1.2

f(x)=x

31

1.3

f(x)=x^2

1.4

f(x)=x^3

32

Bijlage 2: Maple-uitwerking: bepaling coëfficiënten

 23 ℎ 8  3ℎ  23 ℎ = · =0

= ·

= ·

33

Bijlage 3: Maple-uitwerking: bepaling hoogste macht

34

Bijlage 4: Maple-uitwerking: hoogste afgeleide

35

Bijlage 5: Maple-uitwerking: exacte functie-integralen

F:

G:

H:

36

Bijlage 6: Matlab-script: eigenmethode function [benadering, Nf] = eigenmethode(f,a,b,N); %__________________________________________________________________________ ___________________________ % % deze functie inegreerd de ingevoerde functie tussen de (ingevoerde) % grenzen a en b, over N deelintervallen, met behulp van de eigenmethode; % % INVOER % f: de te integreren functie % a: Het beginpunt van het integratie-interval % b: Het eindpunt van het integratie-interval( % N: aantal deelintervallen waarover gintergreerd wordt % % UITVOER % benadering: De benadering van de integratie % Nf: Hierin wordt het aantal functie-evaluaties genoteerd. % % version: 1 date: 1-dec-2010 author: J. Willems %__________________________________________________________________________ ___________________ h = (b-a)/(4*N); k=[0:4*N]; xk = a+k*h; fxk = f(xk); wk =zeros(1,(4*N)+1); wk(1)=2; wk(2:4:4*N)=0; wk(3:4:4*N)=8; wk(4:4:4*N)=0; wk(end)=2; benadering = h/3*(wk*fxk'); Nf = length(xk);

37

Bijlage 7: Matlab-script: Trapeziummethode function [benadering, Nf] = trapeziumregel(f,a,b,N); %__________________________________________________________________________ ___________________________ % % deze functie berekent een benadering van de integraal van de functie f over het interval [ a, b ] % met behulp van de samengestelde trapeziumregel. % INVOER % f: een functie van een variabele, de functie moet kunnen werken op vectoren % a: linker eindpunt van het integratieinterval % b: rechter eindpunt van het integratieinterval % N: het aantal deelintervallen waarin het interval [a,b] wordt verdeeld, en wordt gebruikt om de nauwkeurigheid te sturen % % UITVOER % benadering: in deze variabele wordt de benadering van de integraal afgeleverd % Nf: in deze variabele wordt het aantal functie-evaluatie genoteerd % % voorbeeld van een aanroep: [schatting,aantalF]=trapeziumregel(@(x) x.*sin(3*x.^2), 0, 2, 50) % % % % version: 2.1 date: 13-dec-2007 author: PHM Wolkenfelt %__________________________________________________________________________ ___________________ h = (b-a)/N; k = [0:N]; xk = a+k*h; fxk = f(xk); bijbehorende functiewaarden staan in fxk wk =[1/2 ones(1,N-1) 1/2]; samengestelde trapeziumregel benadering = h*wk*fxk'; is een kolomvector Nf = length(xk);

38

%

xk zijn de basispunten, en de % de gewichten van de % wk is een rij-vector, fxk'

Bijlage 8: Matlab-script: Integraal 1 %-------------------------------------------------------------------------% %

Dit Matlabscript bepaald aan de hand van verschillende numierke methoden de primitieve en integratie van een ingevoerde functie.

% % %

De functie is als volgt gedefinieerd: f(x) = 5*x^4 + 9*x^2 + 4x.

% author: Jeffrey L.A. Willems % version: 2.0 % date: 1-6-2011 %-------------------------------------------------------------------------%%

Rekeneigenschappen in initialisatie van script

clear all clc format long %%

Hoofdberekening van functie

% plot f(x) f=inline( '5.*x.^4 + 9.*x.^2 + 4.*x' );

% functie-implementatie

figure(1) tekengrafiek(f,0,2,0.001,1) title('f(x)=5*(x^4)+9*(x^2)+4*(x)' ) xlabel('x') ylabel('f(x)') % Plot F(x) F=inline( '1.*x.^5 + 3.*x.^3 + 2.*x^2' ); figure(2) tekengrafiek(F,0,2,0.001,2) title('F(x)=1.*x.^5 + 3.*x.^3 + 2.*x^2' ) xlabel('x') ylabel('F(x)') % integraal F(x) van a tot b a=0; b=2; intf=64

%%

% integraal van functie

integratie mbv trapeziummethode

N = [16 32 64 128 256 512 1024] for k=1:7 % de N wordt zo gekozen, dat h=1/8 1/16 1/32 enz.

39

[benaderingTR(k), NfTR(k)] = benadering van de integraal.

trapeziumregel(f,a,b,N(k));

% Exacte fout (1e kolom) [WAREFOUTTR(k,1)] = abs(intf-benaderingTR(k)); verschil tussen de benadering en exacte waarde opgeslagen, dus de ware fout end

wordt

%

% het berekend en

% bepalen van orde en foutschatting for k=1:5; [foutschattingTR(k,1)] = (benaderingTR(k+1)-benaderingTR(k))*(1/3); [ordeTR(k,1)] = log((benaderingTR(k+1)benaderingTR(k))/(benaderingTR(k+2)-benaderingTR(k+1)))/(log(2));

end benaderingTR=benaderingTR' NfTR=NfTR' [WAREFOUTTR] [foutschattingTR] [ordeTR] %% Integratie eigenmethode N=[4 8 16 32 64 128 256]; for k=1:7 % de N wordt zo gekozen, dat h=1/8 1/16 1/32 enz. [benaderingEM(k), NfEM(k)] = benadering van de integraal.

eigenmethode2(f,a,b,N(k));

% Exacte fout (1e kolom) [WAREFOUTEM(k,1)] = abs(intf-benaderingEM(k)); verschil tussen de benadering en exacte waarde opgeslagen, dus de ware fout end

wordt

%

% het berekend en

% bepalen van orde en foutschatting for k=1:5; [foutschattingEM(k,1)] = (benaderingEM(k+1)-benaderingEM(k))*(1/15); [ordeEM(k,1)] = log((benaderingEM(k+1)benaderingEM(k))/(benaderingEM(k+2)-benaderingEM(k+1)))/(log(2));

end benaderingEM=benaderingEM' NfEM=NfEM' [WAREFOUTEM] [foutschattingEM] [ordeEM]

40

%% QUAD en QUADL for t=1:15; tol=1*10^(-t); [BENADERINGQ,function_count]=quad(f,0,2,tol); quadf(t,1)=BENADERINGQ; quadf(t,2)=abs(intf-BENADERINGQ); quadf(t,3)=function_count; [BENADERINGQL,function_count]=quadl(f,0,2,tol); quadlf(t,1)=BENADERINGQL; quadlf(t,2)=abs(intf-BENADERINGQL); quadlf(t,3)=function_count; end quadf quadlf figure(3); loglog(quadf(:,3),quadf(:,2)); xlabel('f-evaluaties' ); ylabel('nauwkeurigheid (ware fout)' ); title('Nauwkeurigheid tegen f-evaluaties' );

figure(4); loglog(quadlf(:,3),quadlf(:,2)); xlabel('f-evaluaties' ); ylabel('nauwkeurigheid (ware fout)' ); title('Nauwkeurigheid tegen f-evaluaties' );

41

% de functie wordt geplot % bijschrift bij x-as % bijschrift bij y-as % titel

% de functie wordt geplot % bijschrift bij x-as % bijschrift bij y-as % titel



% % %

De functie is als volgt gedefinieerd: f(x) = sin((x^4)^(1/7)).





% plot f(x) f=inline( 'sin((x.^4).^(1/7))' );


figure(1) tekengrafiek(f,0,2,0.001,1) title('sin(x^(4/7))' ) xlabel('x') ylabel('f(x)') % integraal F(x) van a tot b a=0; b=1; intf=0.5770813840690 van functie

%%

% integraal


N=[8 16 32 64 128 256 512]; for k=1:7 % de N wordt zo gekozen, dat h=1/8 1/16 1/32 enz. [benaderingTR(k), NfTR(k)] = benadering van de integraal.


% Exacte fout (1e kolom) [WAREFOUTTR(k,1)] = abs(intf-benaderingTR(k)); verschil tussen de benadering en exacte waarde opgeslagen, dus de ware fout end

42

wordt

%

% het berekend en

% bepalen van orde en foutschatting for k=1:5; [foutschattingTR(k,1)] = (benaderingTR(k+1)-benaderingTR(k))*(1/3); [ordeTR(k,1)] = log((benaderingTR(k+1)benaderingTR(k))/(benaderingTR(k+2)-benaderingTR(k+1)))/(log(2));

end benaderingTR=benaderingTR' NfTR=NfTR' [WAREFOUTTR] [foutschattingTR] [ordeTR] %% Integratie eigenmethode N=[2 4 8 16 32 64 128 ]; for k=1:7 % de N wordt zo gekozen, dat h=1/8 1/16 1/32 enz. [benaderingEM(k), NfEM(k)] = benadering van de integraal.



wordt

%

% het berekend en


end benaderingEM=benaderingEM' NfEM=NfEM' [WAREFOUTEM] [foutschattingEM] [ordeEM] %% QUAD en QUADL for t=1:15; tol=1*10^(-t); [BENADERINGQ,function_count]=quad(f,a,b,tol); quadf(t,1)=BENADERINGQ; quadf(t,2)=abs(BENADERINGQ-intf); quadf(t,3)=function_count;

43

[BENADERINGQL,function_count]=quadl(f,a,b,tol); quadlf(t,1)=BENADERINGQL; quadlf(t,2)=abs(BENADERINGQL-intf); quadlf(t,3)=function_count; end quadf quadlf figure(3); loglog(quadf(:,3),quadf(:,2),quadlf(:,3),quadlf(:,2)); de functie wordt geplot xlabel('f-evaluaties' ); % bijschrift bij x-as ylabel('nauwkeurigheid (ware fout)' ); % bijschrift bij y-as title('Nauwkeurigheid tegen f-evaluaties QUAD en QUADL' ); % titel

44

%



% % %

De functie is als volgt gedefinieerd: f(x) = sin((x^4)^(1/7)).





% plot f(x) f=inline( 'sin((x.^4).^(1/7))' );


figure(1) tekengrafiek(f,0,2,0.001,1) title('sin(x^(4/7))' ) xlabel('x') ylabel('f(x)') % integraal F(x) van a tot b a=1; b=2; intf=0.9412826152700 %%

% integraal van functie


N=[8 16 32 64 128 256 512]; for k=1:7 % de N wordt zo gekozen, dat h=1/8 1/16 1/32 enz. [benaderingTR(k), NfTR(k)] = benadering van de integraal.


% Exacte fout (1e kolom) [WAREFOUTTR(k,1)] = abs(intf-benaderingTR(k)); verschil tussen de benadering en exacte waarde opgeslagen, dus de ware fout end % bepalen van orde en foutschatting for k=1:5;

45

wordt

%

% het berekend en

[foutschattingTR(k,1)] = (benaderingTR(k+1)-benaderingTR(k))*(1/3); [ordeTR(k,1)] = log((benaderingTR(k+1)benaderingTR(k))/(benaderingTR(k+2)-benaderingTR(k+1)))/(log(2)); end benaderingTR=benaderingTR' NfTR=NfTR' [WAREFOUTTR] [foutschattingTR] [ordeTR] %% Integratie eigenmethode N=[2 4 8 16 32 64 128]; for k=1:7 % de N wordt zo gekozen, dat h=1/8 1/16 1/32 enz. [benaderingEM(k), NfEM(k)] = benadering van de integraal.



wordt

%

% het berekend en


end benaderingEM=benaderingEM' NfEM=NfEM' [WAREFOUTEM] [foutschattingEM] [ordeEM] %% QUAD en QUADL for t=1:15; tol=1*10^(-t); [BENADERINGQ,function_count]=quad(f,a,b,tol); quadf(t,1)=BENADERINGQ; quadf(t,2)=abs(BENADERINGQ-intf); quadf(t,3)=function_count; [BENADERINGQL,function_count]=quadl(f,a,b,tol); quadlf(t,1)=BENADERINGQL; quadlf(t,2)=abs(BENADERINGQL-intf);

46

Opdracht Numerieke Methoden; Numerieke Integratie Juni 2011

Recommend Documents