Numerieke wiskunde: Oefeningen

Numerieke wiskunde: Oefeningen Philippe Nimmegeers 2de bachelor ingenieurswetenschappen Katholieke Universiteit Leuven Academiejaar: 2010-2011

Numerieke Wiskunde: Oefeningen Philippe Nimmegeers

1

Inhoudsopgave Oefenzitting 1: Foutenanalyse................................................................................................................ 4 Oefening 1: Conversie van binair naar decimaal................................................................................. 4 Oefening 2: Conversie van decimaal naar binair................................................................................. 4 Oefening 3: Voorstelling van getallen in Matlab ................................................................................. 5 Oefening 4: Berekenen van

................................................................................. 6

Oefening 5: Berekenen van een machtreeks ...................................................................................... 9 Extra 1: Een onschadelijke berekening? ............................................................................................ 10 Extra 2: Berekening van

................................................................................................. 15

Extra 3: Afrondingsfouten zijn niet random (W. Kahan) ................................................................... 16 Extra 4: Een benadering van

.......................................................................................................... 17

Extra 5: Het bepalen van de basis b en het aantal cijfers in de mantisse p ...................................... 18 Oefenzitting 2: Het oplossen van stelsels lineaire vergelijkingen ......................................................... 19 Oefening 1: De LU-ontbinding ........................................................................................................... 19 Oefening 2: De LU-ontbinding, deel twee ......................................................................................... 22 Oefening 3: Conditie en stabiliteit..................................................................................................... 23 Oefening 4: Een meetkundige interpretatie van de conditie van een stelsel in 2D.......................... 25 Oefening 5: Rij-schalering.................................................................................................................. 26 Oefenzitting 3: Veelterminterpolatie .................................................................................................... 27 Oefening 1: Interpolatie volgens Lagrange ....................................................................................... 27 Oefening 2: Interpolatie volgens Newton ......................................................................................... 31 Oefening 3: Conditie van het interpolatieprobleem ......................................................................... 33 Oefening 4: Multiple choice .............................................................................................................. 34 Oefenzitting 4: Numerieke differentiatie en integratie ........................................................................ 35 Oefening 1: Numerieke differentiatie ............................................................................................... 35 Oefening 2: Numerieke integratie: Enkelvoudige interpolerende kwadratuurformule ................... 37 Oefening 3: Automatische integratie met de middelpuntsregel ...................................................... 38 Oefenzitting 5: Differentiaalvergelijkingen ........................................................................................... 39 Oefening 1: De slinger ....................................................................................................................... 39 Oefening 2: Orde van een methode bepalen .................................................................................... 42 Oefening 3: Stabiliteitsgebieden ....................................................................................................... 43 Oefenzitting 6: Het iteratief oplossen van een niet-lineaire vergelijking ............................................. 44 Oefening 1: Numerieke experimenten met bissectie, secantmethode, Newton-Raphson-methode ........................................................................................................................................................... 44


2

Oefening 2: Consistentie en convergentie van stationaire een-staps-substitutiemethoden ........... 49 Oefening 3: Waar of niet waar .......................................................................................................... 54 Oefening 4: Newton-Raphson methodes voor het berekenen van ............................................... 55 Extra 1: Regula falsi ........................................................................................................................... 61 Oefenzitting 7: Het iteratief oplossen van stelsels (niet-) lineaire vergelijkingen ................................ 62 Oefening 1: Conditie van een wortel van een niet-lineaire vergelijking ........................................... 62 Oefening 2: Stelsels niet-lineaire vergelijkingen ............................................................................... 63 Oefening 3: Convergentie voor de methode van Gauss-Seidel voor stelsels lineaire vergelijkingen65 Oefenzitting 8: Eigenwaarden, Optimalisatie en Herhaling ................................................................. 66 Oefening 1: Het berekenen van eigenwaarden ................................................................................ 66 Oefening 2: Optimalisatie: waar of niet waar? ................................................................................. 67 Oefening 3: Herhaling: Stelsels lineaire vergelijkingen .................................................................... 68 Oefening 4: Herhaling: Stelsels lineaire vergelijkingen: convergentie van iteratieve methoden ..... 69 Oefening 5: Herhaling: Stabiliteit van de methode van Euler ........................................................... 71 Examenvraag 1: LU-factorisatie......................................................................................................... 72 Examenvraag 2: Juiste cijfers ............................................................................................................ 74 Examenvraag 3: waar of niet waar? .................................................................................................. 76 Examenvraag 4: Stabiliteit van methoden voor het oplossen van differentiaalvergelijkingen ........ 77 Examenvraag 5: Stelsels differentiaalvergelijkingen (Nog te maken) ............................................... 79 Examenvraag 6: Optimalisatie in meerdere veranderlijken ............................................................. 82 Oefenzitting 9: Optimalisatie in meerdere veranderlijken .................................................................. 85 Oefening 1: De methode van de steilste afdaling en methode van Newton .................................... 85 Oefening 2: Vergelijking methode van Newton en steilste afdaling ................................................. 87 Oefening 3: Methode der steilste afdaling ....................................................................................... 88 Oefening 4: Optimalisatie van een functie in twee veranderlijken. ................................................. 91


3

Oefenzitting 1: Foutenanalyse Oefening 1: Conversie van binair naar decimaal Er wordt gebruik gemaakt van de klassieke getallenvoorstelling met grondtal of basis 2, die door onderstaande uitdrukking wordt weergegeven: ∑

∑

Dit levert onderstaande resultaten op voor de conversie van de gegeven binaire getallen, naar het decimaal talstelsel.

Oefening 2: Conversie van decimaal naar binair Zoals in de opgave aangegeven wordt er gebruik gemaakt van onderstaande website om de voorstelling van een bepaald getal in verschillende talstelsels te geven.


4

Oefening 3: Voorstelling van getallen in Matlab In de cursus vindt men machinenauwkeurigheid:

onderstaande

formule

voor

het

berekenen

van

de

Bereken de machinenauwkeurigheid voor Matlab, indien gegeven is dat men werkt met een basis gelijk aan twee en p gelijk is aan 53. Vergelijk deze waarde met de waarde van eps(1). 

Wat is het verband tussen de machinenauwkeurigheid en eps(1)?

De machinenauwkeurigheid is een bovengrens voor de relatieve fout die we maken als we een getal x door fl(x) vervangen. Voor eps(1) kan men lezen dat deze wordt gedefinieerd als de afstand tussen het getal 1 en het eerstvolgende (grotere)machinegetal. De absolute fout die je maakt door een getal x, met 1 < x en x < 1+ eps(1), door fl(x) voor te stellen is het grootst in het midden van het interval [1, 1 + eps(1)]. Rond 1 is de absolute fout nagenoeg gelijk aan de relatieve fout dus geldt dat:



Hoeveel is (1+0.60 eps)-1? Verklaar de uitkomst!

In Matlab is eps gelijk aan eps(1). Alle getallen die tussen 1 en 1+eps(1) liggen, kunnen niet exact voorgesteld worden. Bijgevolg wordt het getal afgerond naar het getal . Het resultaat van bovenstaande uitdrukking kan dus geschreven worden als: (


)

5

Oefening 4: Berekenen van

(

)

Stel dat we numeriek een limiet willen berekenen: (

)

Om de grafiek te potten van de absolute fout op een semilogplot voor opeenvolgende machten van 10, kan men volgend algoritme schrijven in Matlab. % % % % %

Deze methode tekent de grafiek van de absolute fout op een semilogplot. lim (x--> +oneindig) (1+1/x)^x = exp(1) Bovenstaande uitdrukking wordt getest met opeenvolgende machten (1ste macht t.e.m. 20ste macht) van 10

for i = 1:20 x = 10.0^i; absolutefout(i) = (1+1/x)^x; end; % plotten van de grafiek van de absolute fout figure(1) semilogy(abs(absolutefout-exp(1))) title('absolute fout')

Dit levert onderstaande foutencurve op:


6

 Bij

Waarom blijft de fout constant vanaf 1016? of groter telt Matlab een getal op bij 1 dat kleiner is dan de

machinenauwkeurigheid (

). Bijgevolg zal Matlab dit getal afronden naar 1, zodat

de fout constant blijft. 

Waarom stijgt de fout vanaf 108?

Vanaf zijn de afrondingsfouten groter dan de benaderingsfouten. Dit kan men zien op de foutcurve. Bij zijn er geen fouten en zijn er nog 16 cijfers juist. Bij enkel cijfer meer juist.

is er geen

In het begin van de grafiek is enkel de benaderingsfout van belang. Op de grafiek kan men zien dat de benaderingsfout lineair afneemt (op logaritmische schaal!) Aan het tweede deel van de grafiek (waar dus de afrondingsfouten groter zijn dan de benaderingsfouten) kan een lineaire toename waargenomen worden. Het afnemen van het aantal juiste cijfers en dus de afrondingsfouten hebben een lineair verloop (op logaritmische schaal!). De absolute fout is de som van de benaderingsfout en de afrondingsfout en geeft bovenstaande grafiek weer. Deze lijnen worden gecombineerd en snijden elkaar in het punt waar , wat ervoor zorgt dat de fout stijgt vanaf deze waarde voor n. Men zou kunnen proberen om het resultaat te verbeteren door de limiet op een andere manier uit te rekenen: ( 

(

))

Waarom haalt het berekenen van de limiet op bovenstaande manier niets uit? Er is een klein verschil tussen de manieren waarop de limiet berekend wordt. Er zijn nl. andere afrondingsfouten. Het algemeen verloop van de foutencurve is echter analoog. Het berekenen van de logaritme en de exponentiele leveren ook afrondingsfouten op. Het probleem is echter dat de bewerking 1+1/n steeds wordt uitgevoerd.

Deze methode kan verbeterd worden door de reeksontwikkeling in te vullen voor

(

): ( 

(

))

( [

])

Wat is het probleem met deze formule?


7

Aangezien men twee grote getallen, die ongeveer even groot zijn, van elkaar gaat aftrekken, wordt de relatieve fout enorm opgeblazen. Blijkbaar hebben we in basis 2 een veel beter resultaat dan in basis 10. 

Welk probleem hebben we nu gedeeltelijk vermeden? Men verkrijgt door het nemen van basis 2 minder afrondingsfouten, aangezien de machten van 2 voorstelbaar zijn, zolang deze binnen het bereik van de bewegende kommavoorstelling blijven.


8

Oefening 5: Berekenen van een machtreeks De machtreeks van exp(x) rond x=0, wordt gegeven door: ∑ Deze machtreeks kan worden geëvalueerd voor x=1 en afgebroken na een aantal termen en dit geeft dan een benadering voor het getal e. 

Hoeveel termen zijn nodig om op deze manier de beste benadering te bekomen? Waarom draagt volgende term niet meer bij aan het geheel? Er zijn 8 termen nodig (fout[7] is minimaal). De volgende term is gelijk aan , wat kleiner is dan . Bijgevolg levert de volgende term geen bijdrage meer, aangezien er telkens “0” wordt bijgeteld.



Waarom wordt de fout met vier termen fout[3] groter als je met een precisie van zes cijfers werkt? Er zijn afrondingsfouten ten gevolge van er meer digits worden beschouwd.



Geef een mogelijkheid om de som nauwkeuriger uit te rekenen in een gegeven precisie. Men kan de termen optellen van klein naar groot, wat zorgt voor meer nauwkeurigheid.


9

Extra 1: Een onschadelijke berekening? Methode om de matrix A op te stellen die gegenereerd wordt door de opgegeven lus te laten doorvoeren. function [A] = lusje1(x) %UNTITLED Summary of this function goes here % Detailed explanation goes here for i = 1 : 40 x = sqrt(x); A(i,1) = x; end for i = 1 : 40 x = x^2; end A end

Methode om de matrix verschil te berekenen. function [D] = bewerking1advanced(A,B) %UNTITLED3 Summary of this function goes here % Detailed explanation goes here C = A-B; for i = 1 : 40 x = abs(C(i,1)); D(i,1) = x; end end A = lusje1(100); verschil = bewerking1advanced(A, exacte_x); plot(verschil,’b’)


10

Omgekeerde lus: kwadrateren function [A] = omgekeerdlusje1(x) %UNTITLED4 Summary of this function goes here % Detailed explanation goes here for i = 1 : 40 x = x^2; A(i,1) = x; end for i = 1 : 40 x = sqrt(x); end end

Vervolgens worden volgende commando’s aangeroepen: A = omgekeerdlusje1(100); verschil = berekening1advanced(A, exacte_x;) semilogy(verschil,'b')


11

Voor i van 1 tot en met 60: function [A] = lusje2(x) %UNTITLED5 Summary of this function goes here % Detailed explanation goes here for i = 1 : 60 x = sqrt(x); A(i,1) = x; end for i = 1 : 60 x = x^2; end end function [D] = bewerking2advanced(A,B) %UNTITLED3 Summary of this function goes here % Detailed explanation goes here C = A-B; for i = 1 : 60 x = abs(C(i,1)); D(i,1) = x; end end exacte = [exacte_x; ones(20,1)]; A = lusje2(100); verschil = bewerking2advanced(A, exacte); plot(verschil,’b’)

exacte = [exacte_x; ones(20,1)]; A = omgekeerdlusje2(100); verschil = bewerking2advanced(A, exacte); semilogy(verschil,’b’)


12

Omgekeerde lus: kwadrateren function [A] = omgekeerdlusje2(x) %UNTITLED4 Summary of this function goes here % Detailed explanation goes here for i = 1 : 60 x = x^2; A(i,1) = x; end for i = 1 : 60 x = sqrt(x); end end exacte = [exacte_x; ones(20,1)]; A = omgekeerdlusje2(100); verschil = bewerking2advanced(A, exacte); semilogy(verschil,’b’)


13

Aangezien bovenstaande functions gekend zijn, kunnen onderstaande functies gedefinieerd worden en kan er eenvoudig voor de verschillende gevallen nagegaan worden hoe de grafieken eruit zien. function [] = Vierkantswortel(x) %UNTITLED7 Summary of this function goes here % Detailed explanation goes here load exacte_x.mat exacte = [exacte_x; ones(20,1)]; A = lusje2(x); verschil = bewerking2advanced(A,exacte); if x<10 plot(verschil,'b') else semilogy(verschil,'b') end end function [] = kwadrateren(x) %UNTITLED6 Summary of this function goes here % Detailed explanation goes here load exacte_x.mat exacte = [exacte_x; ones(20,1)]; A = omgekeerdlusje2(x); verschil = bewerking2advanced(A,exacte); if x<10 plot(verschil,'b') else semilogy(verschil,'b') end end

Vb: x = 0.5 Vierkantswortel(0.5)

kwadrateren(0.5)


14

Extra 2: Berekening van function [f] = exp_1(x) %UNTITLED Summary of this function goes here % Detailed explanation goes here if x==0 f = 1; else f = (exp(x)-1)/x; end end function [f] = exp_2(x) %UNTITLED2 Summary of this function goes here % Detailed explanation goes here y = exp(x); if y==1 f = 1; else f = (y-1)/log(y); end end


15

Extra 3: Afrondingsfouten zijn niet random (W. Kahan) function r=not_rand %Uit Nicholas J. Higham "Accuracy and Stability of Numerical Algorithms" for i=1:300 x=1.606+(i-1)*2^(-52); r(i)=horner(x,[4,-59,324,-751,622])/horner(x,[1,-14,72,-151,112]); end plot(r-8.75237658077849,'.'); %% function y=horner(x,p); n=length(p); y=p(1); for i=2:n y=y*x+p(i); end


16

Extra 4: Een benadering van 

10000 termen



Cfr 1a.


17

Extra 5: Het bepalen van de basis b en het aantal cijfers in de mantisse p function b = bepaal_b %BEPAAL_B bepaalt het grondtal waarmee de computer werkt % % B = BEPAAL_B retourneert het grondtal waarmee de computer werkt. a = 1; while ((a + 1) - a) == 1 a = a * 2; end; i = 1; while a == (a + i) i = i + 1; end; b = (a + i) - a; end function p = bepaal_p %BEPAAL_P bepaalt het aantal beduidende cijfers in de mantisse. % % P = BEPAAL_P retourneert het aantal beduidende cijfers in de mantisse.

% Voor dit algoritme hebben we de basis nodig, % we gebruiken hiervoor het bepaal_b script. b = bepaal_b; p = 1; z = b; while ((z + 1) - z) == 1 p = p + 1; z = z * b; end end

De basis is gelijk aan 2 en p gelijk aan 53


18

Oefenzitting 2: Het oplossen van stelsels lineaire vergelijkingen Oefening 1: De LU-ontbinding De LU-ontbinding van een matrix gebeurt in feite impliciet wanneer je een stelsel oplost aan de hand van Gauss-eliminatie. 1. Stel dat volgende matrix gegeven is: [ 

]

Wat is LU? L is een eenheidsbenedendriehoeksmatrix die de inverse is van de benedendriehoeksmatrix T die indien links vermenigvuldigd met A een bovendriehoeksmatrix U oplevert. Bijgevolg is LU gelijk aan de matrix A.



Wat zijn de Matlab-bevelen om uit L en U de determinant van A te berekenen? De determinant van A kan eenvoudig berekend worden door de determinant van LU te berekenen. Aangezien de determinant van een product gelijk is aan het product van de determinanten, kan de determinant van A geschreven worden d.m.v. onderstaande uitdrukking:

Aangezien men weet dat L een eenheidsbenedendriehoesmatrix is en de determinant van een benedendriehoeksmatrix gelijk is aan het product van de diagonaalelementen, is de determinant van L gelijk aan 1. Bijgevolg kan de determinant van A geschreven worden als:

Aangezien de matrix U een bovendriehoeksmatrix is, is de determinant van U gelijk aan het product van de diagonaalelementen en kan men eenvoudigweg in Matlab, in plaats van de determinant zelf van U te berekenen, het product van de diagonaalelementen van U berekenen. Dit levert onderstaande Matlab-code op voor het berekenen van de determinant van A : (

)

Er dient opgemerkt te worden dat bovenstaande uitwerking enkel geldig is voor gevallen waarbij er geen rijverwisselingen worden doorgevoerd. Indien men te maken heeft met rijpivotering dient de permutatiematrix P nog in rekening gebracht te worden.


19

In dit geval is de determinant van A op het teken na bepaald. De determinant van een permutatiematrix kan gelijk zijn aan +1 en de determinant van de eenheidsbenedendriehoeksmatrix L blijft gelijk aan 1. 

Stel dat een rechterlid b gegeven is. Hoe kan men dan m.b.v. de matrices L en U het stelsel Ax = b eenvoudig oplossen? Aangezien LU gelijk is aan A, kan het gegeven stelsel ook geschreven worden als:

Daar een matrixproduct uitrekenen even veel rekentijd in beslag neemt als het berekenen van een LU-ontbinding, dient men voorzichtig te werk te gaan. Men kan eerst stellen dat Ux = y, zodat bovenstaand stelsel geschreven kan worden als volgt: y=forward(L,b)

Bovenstaand benedendriehoekig stelsel kan eenvoudig opgelost worden naar y. Aangezien men y bekomen heeft, kan x gevonden worden door volgend stelsel op te lossen. x=backsubs(U,y)

Dit komt eveneens neer op het oplossen van een eenvoudig stelsel, dat in dit geval wel een bovendriehoekig stelsel is. 2. Stel dat onderstaande matrix gegeven is: [ 

]

Waarom is L geen benedendriehoeksmatrix meer? Het spilelement a33=0=u33. Er worden bijgevolg bij het zoeken van de beste pivot rijverwisselingen doorgevoerd. De rijverwisselingen zijn verwerkt in de matrix L die het algoritme doorvoert.



Bekijk de Matlab-documentatie van lu nog eens en gebruik de vorm [L, U, P] = lu(A). Wat is LU nu? PA=LU


20



Hoe ga je in dit geval het stelsel Ax=b oplossen? Aangezien Ax=b, kunnen beide leden vermenigvuldigd worden met de matrix P. bijgevolg bekomt men een nieuw stelsel:

Aangezien LU= PA, kan dit stelsel herschreven worden als volgt:

Stel dat Ux =y, dan bekomt men onderstaand benedendriehoekig stelsel, dat kan opgelost worden met voorwaartse substitutie:

Aangezien de matrix y, bepaald wordt uit bovenstaand stelsel, dient x bepaald te worden door het oplossen van onderstaand bovendriehoekig stelsel:


21

Oefening 2: De LU-ontbinding, deel twee Genereer een stelsel van zes vergelijkingen met het script genstelsel.m. Je bekomt een A, een b en de exacte oplossing x. Met behulp van de functies lutx.m en lutx2.m en daarna de functies forward.m en backsubs.m kan dit stelsel opgelost worden. 

Vergelijk nu de relatieve fout van de twee bekomen oplossingen. Is het verschil te verklaren door de stabiliteit van het algoritme of door de conditie van het probleem? Hetzelfde probleem wordt opgelost met twee verschillende algoritmes en de relatieve fouten verschillen in grootteorde in beide gevallen. De verschillende grootteorde van de fout is bijgevolg te verklaren door de stabiliteit van de algoritmes.



Wat doet lutx anders dan lutx2? lutx gebruikt het eerste niet-nulelement als pivot, terwijl lutx2 op zoek gaat naar het grootste niet-nulelement.


22

Oefening 3: Conditie en stabiliteit We genereren nu opnieuw een stelsel van zes vergelijkingen, maar ditmaal met het script genstelselc.m, je bekomt opnieuw een matrix A, een rechterlid b en de exacte oplossing x. Los nu opnieuw dit stelsel op met behulp van lutx en lutx2 en nadien forward en backsubs. Noem de oplossingen opnieuw y1 en y2. 

Vergelijk nogmaals de relatieve fout van de twee bekomen oplossingen. Wat merk je op? Ligt dit aan de stabiliteit van het algoritme of aan de conditie van het probleem? De twee methodes leveren eenzelfde even grote fout op voor hetzelfde probleem. De relatieve fout is het gevolg van de slechte conditie van het probleem.



Wat is het conditiegetal van de matrix A?



‖ ‖‖ ‖ In de cursus vind je volgende uitdrukking

‖

‖

‖ ‖

‖

‖

‖ ‖

terug m.b.t het

conditiegetal van de matrix A. Kan je aan de hand hiervan de relatieve fouten voorspellen? ( ) Voor

‖

‖

‖ ‖

kan men stellen dat dit ongeveer gelijk is aan de machinenauwkeurigheid,

aangezien A tot op de machinenauwkeurigheid gekend is. Dit levert op, mits invullen in de formule:



Genereer nu twee stelsels van 5 vergelijkingen met de volgende commando's: [M1, x1] = genstelseld(5,0) en [M2, x2] = genstelseld(5, 1e-10). Bereken opnieuw de relatieve fouten en vul verder aan:



Verklaar de verschillen. Heeft dit te maken met de stabiliteit van het algoritme of de conditie van het probleem? De problemen zijn goed geconditioneerd. De verschillen zijn te verklaren door de stabiliteit van de algoritmes.


23



Wat is de functie van de parameter delta (die je voor M1 op 0 hebt gezet en voor M2 op 1010)? (Test dit eens met pen en papier op een 3x3 stelsel.)


24

Oefening 4: Een meetkundige interpretatie van de conditie van een stelsel in 2D Op pagina 65 van de cursus wordt voor een 2x2 stelsel de conditie van een stelsel grafisch voorgesteld. We zullen nu in Matlab experimenteren met twee stelsels. Gebruik hiervoor het script tekenstelsel.m. Met de parameters alpha en beta kan je een perturbatie aanbrengen op de matrix A (althans de eerste kolom van A) en op het rechterlid b. Deze perturbaties worden als volgt aangebracht: en .



Wat is de conditie van beide matrices A? Is dit goed of slecht? Probeer ook eens cond(eye(2)). Bereken ook eens det (A). Beide stelsels zijn bijgevolg goed geconditioneerd. cond(eye(2))=1 det(A1)=2, det(A2)=5.


25

Oefening 5: Rij-schalering Pivotering helpt niet om slecht geconditioneerde problemen op te lossen. Slecht geconditioneerde problemen kan men wel herformuleren, zodat de conditie verbetert. Bijvoorbeeld door het probleem of een deel van het probleem te herschalen. Stel we hebben volgend orthogonaal steslel:



Voor welke twee mogelijke problemen waarschuwt Matlab?



Kan je de twee mogelijke oorzaken uitleggen? Close to singular betekent dat de rijen van de matrix bijna lineair afhankelijk zijn van elkaar. (~singuliere matrix) Badly scaled betekent dan weer dat de matrixelementen van heel verschillende grootteordes zijn.



Wat kan je doen met dit stelsel om de conditie te verbeteren in het geval van slechte schalering? Door de tweede vergelijking te delen door , bekomt men een stelsel met beter geschaalde matrix en bijgevolg ook een beter geconditioneerd probleem.



Wat is de conditie van het verbeterde probleem? De conditie van het gewijzigde stelsel is gelijk aan 1.



Kun je in het andere geval (bijna singulier) de conditie verbeteren? Nee.


26

Oefenzitting 3: Veelterminterpolatie Oefening 1: Interpolatie volgens Lagrange Gegeven een tabel met koppels geschreven worden als:

kan de interpolerende veelterm van graad n

De veeltermen li(x) worden Lagrangeveeltermen genoemd. De bijhorende formule voor yn(x) is de interpolerende veelterm volgens Lagrange. 1. Bewijs onderstaande uitdrukkingen: 

Bewijs dat :

Bewijs: Stel Dan:

∑

De interpolatiefout is gelijk aan het verschil van de functie zelf en de interpolerende veelterm en kan ook beschreven worden door onderstaande uitdrukking:

Bovenstaande uitdrukking is gelijk aan nul, aangezien k n. Bijgevolg is de interpolatiefout gelijk aan nul en geldt het volgende:

Hieruit volgt dat:

Dit is de uitdrukking die bewezen moest worden.


27



Bewijs dat:

Bewijs: Stel Dan:

∑

De interpolatiefout is gelijk aan het verschil van de functie zelf en de interpolerende veelterm en kan ook beschreven worden door onderstaande uitdrukking:

Bovenstaande uitdrukking is gelijk aan nul. Bijgevolg is de interpolatiefout gelijk aan nul en geldt het volgende:

Hieruit volgt dat:

Dit is de uitdrukking die bewezen moest worden.


28

2. Gegeven de equidistante interpolatiepunten:

Gebruik gegeven door:



als nieuwe veranderlijke. Dan worden de Lagrangeveeltermen

Wat is het verband tussen u en k? Men kan k schrijven in functie van u door de gegeven uitdrukkingen als volgt in elkaar te substitueren: (

)

(

)

In de interpolatiepunten volgt uit bovenstaande uitdrukking, dat k=u, aangezien x de waarde xk kan aannemen.  Volgens de gegeven formule voor de Lagrangeveeltermen geschreven worden als:

Lagrangeveeltermen,

kunnen

de

̌ ̌ ̌ De Lagrangeveelterm ̌ onderstaande resultaten op:

kan geëvalueerd in de opgegeven punten. Dit levert ̌ ̌ ̌ ̌ ̌

Uit bovenstaande resultaten kan men vaststellen dat in alle interpolatiepunten het volgende geldt: ̌


29

3. Benader onderstaande derdegraadsveelterm door een tweedegraadsveelterm op [0,3], neem als interpolatiepunten 0,1 en 3. Maak ook een schatting voor de interpolatiefout En(x) in dit interval (n= 2,3). Eerst wordt de derdegraadsveelterm geëvalueerd in de interpolatiepunten. Dit levert onderstaande functiewaarden op:

De tweedegraadsveelterm kan opgesteld worden d.m.v. onderstaande uitdrukking:

Eerst worden

en

berekend.

Vervolgens worden de Lagrangeveeltermen geschreven.

De interpolerende tweedegraadsveelterm kan bijgevolg geschreven worden als:

De interpolatiefout van de tweedegraadsveelterm kan geschreven worden d.m.v. de definitie:

Indien men een derdegraadsveelterm zal beschouwen als interpolerende veelterm, zal men En=0 uitkomen voor de interpolatiefout, aangezien de derdegraadsveelterm benaderd wordt door een derdegraadsveelterm kan er immers geen fout gemaakt worden.


30

Oefening 2: Interpolatie volgens Newton 1. Gegeven een tabel van meetwaarden bekomen door de functie bemonsteren (samplen) in de x-waarden: 0, 1, 2, 3 en 4.

te

Stel je wil een veelterminterpolatie berekenen volgens Newton. Schrijf het stelsel Ax=b neer dat je dan moet oplossen om de coëfficiënten te bepalen. Een veelterminterpolatie volgens Newton kan berekend worden als volgt:

Dit levert onderstaand stelsel op:

[

𝑥 𝑥 𝑥 𝑥

𝑥 𝑥 𝑥 𝑥

𝑥 𝑥 𝑥

𝑥 𝑥 𝑥

𝑥 𝑥 𝑥

𝑥 𝑥 𝑥

𝑥 𝑥

𝑥 𝑥

𝑥 𝑥

𝑥 𝑥

𝑥 𝑥

𝑥 𝑥

𝑥

𝑥

𝑥

𝑥

𝑥

𝑥

𝑥

𝑏 𝑏 𝑏 𝑏 𝑥 ] [𝑏 ]

[

]

Dit stelsel kan met ingevulde waarden voor xi geschreven worden als:

[


][ ]

[

]

31

2. Een veelterm kan in verschillende basissen voorgesteld worden:

Bekijk de grafieken van deze basissen voor de veeltermen van graad 5. Geef een intuïtieve verlaring voor de uitspraak “De coëfficiënten van de veelterm ten opzichte van de monomiale basis zijn (voor hoge graad n) gevoelig voor de functiewaarden f i (d.w.z. de bepaling van de coëfficiënten is een slecht gesteld probleem indien de graad n hoog is). Voor toenemende graad wordt het steeds moeilijker om de functies van elkaar te onderscheiden, wat ervoor zorgt dat de kolommen van de Vandermonde-matrix bijna lineair afhankelijk worden. Voor de meeste keuzes van xi groeit het conditiegetal van een Vandermonde-matrix minstens exponentieel met het aantal gegevens n+1.


32

Oefening 3: Conditie van het interpolatieprobleem Een monomiale basis wordt voorgesteld door onderstaande matrix (Vandermondematrix):

Een Newtonbasis kan opgesteld worden door onderstaande benedendriehoeksmatrix:

De Lagrangebasis kan men voorstellen d.m.v. een eenheidsmatrix. Bij een stijgende graad zijn de matrices van beide methoden slechter geconditioneerd, maar de matrix van de Newton methode is steeds net iets beter. De Lagrangebasis wordt voorgesteld door een eenheidsmatrix, waarvan men weet dat het conditiegetal gelijk is aan 1. Bij punten zeer dicht bij elkaar gaat de determinant van de Vandermondematrix naar nul. ∏


33

Oefening 4: Multiple choice

1. Waar 2. Niet waar: x0(0,0), x1(1,2) =>graad 0: er is geen enkele constante die door beste waarde gaat. 3. Waar: zie stelling (ze bestaat en is enig) 4. Waar: Alle veeltermen van graad n-1 en hoger. 5. Niet waar, vb.: boogtangensfunctie.


34

Oefenzitting 4: Numerieke differentiatie en integratie Oefening 1: Numerieke differentiatie 1. Bewijs onderstaande formules:

Oplossing: 

(Tweede orde Taylorveelterm van x+h) (Veranderen van lid) (Bovenstaande uitdrukking herschreven naar f’(x)) (

)



(Tweede orde Taylorveelterm van x-h) (Veranderen van lid) (Bovenstaande uitdrukking herschreven naar f’(x)) (

)



(Derde orde Taylorveelterm x+h) (“ “ “ x-h) (Aftrekken Taylorveeltermen) (Delen door 2h) (Herschrijven naar f’(x)) (



Analoog als hierboven: echter vierde orde Taylorveeltermen opstellen, herschrijven naar f”(x)/2 en optellen.


35

2. Volgens bovenstaande formules wordt de benaderingsfout kleiner als h kleiner wordt. Toch mag h niet zeer klein gekozen worden. Voor formule 1 en 2 bijvoorbeeld, wordt aangeraden om , met √ de machinenauwkeurigheid, te kiezen. √ Waarom? Hint: je mag veronderstellen dat f(x) en f(x+h) van orde 1 (O(1)) zijn. Oplossing: Een computer of rekenmachine berekent d;m.v. floating point benaderingen (met : (

)

(

[

)

(

) ]

(

)

(

)

3. Stel dat je zelfs het functievoorschrift van f(x) niet kent, maar enkel een aantal functiewaarden in tabelvorm. Hoe kan je dan de afgeleide(n) van f(x) in x=xi numeriek berekenen? Mogen de {xi}ni:0 volledig willekeurig gekozen worden? De afgeleiden kunnen berekend worden door Taylorveeltermbenaderingen te schrijven, zoals aangegeven in puntje 1. Het is echter wel van belang dat men werkt met equidistante punten, opdat men de afgeleide van f(x) zou kunnen bepalen.


36

Oefening 2: Numerieke integratie: Enkelvoudige interpolerende kwadratuurformule Gegeven de niet-equidistante abscissen waarbij h zodanig gekozen is dat enkelvoudige interpolerende kwadratuurformule op om de integraal ∫

Stel een te benaderen:

Eerst wordt een kwadratuurformule opgesteld d.m.v. de interpolerende veeltermen van Lagrange. De gewichten worden in dit geval gegeven door volgende uitdrukking: ∫ Bijgevolg kunnen de lagrangeveeltermen opgesteld worden:

Voor de methode van de onbepaalde coëfficiënten bestaat uit het oplossen van onderstaand stelsel:

{

Aangezien de momenten gegeven worden door onderstaande uitdrukking: ∫

bekomt men onderstaand stelsel:

{


37

Oefening 3: Automatische integratie met de middelpuntsregel Als het aantal deelintervallen in elke stap verdrievoudigt, wordt de pas voor de k-de stap gegeven door:

Bij de middelpuntsregel ligt de eerste abscis De abscissen voor de (k+1)-ste stap kunnen dan geschreven worden als:

Hieruit volgt dat

En de (k+1)-ste benadering kan als volgt ontbonden worden: ( (

)

∑( ( ∑ (

(

) )

∑( (

(

)

∑( ( )

))

(


( )

)) (

))

))

38

Oefenzitting 5: Differentiaalvergelijkingen Oefening 1: De slinger function rechterlid = slinger(t,theta_en_omega) % SLINGER(t,theta_en_omega) = geeft het rechterlid van % de differentiaalvergelijking van de slinger % % Invoer: % t = de tijdsvariabele % theta_en_omega = een vector met 2 elementen : % theta_en_omega(1) = theta % theta_en_omega(2) = omega % % Uitvoer: % rechterlid = kolomvector (!) met 2 elementen % % Info: % De differentiaalvergelijking % theta'' + (c/(m*l))*theta' + (g/l)*sin(theta) = 0 % die de slinger beschrijft wordt herschreven in de vorm van % een stelsel van 2 differentiaalvergelijkingen door het % invoeren van omega := theta'. Na substitutie van a := g/l % en b := c/(m*l) (zie opgaveblad) bekomen we : % { theta' = omega % { omega' = -a*sin(theta) - b*omega % % Deze matlabfunctie - die je zelf moet aanvullen % geeft het rechterlid van deze differentiaalvergelijking % als theta en omega (en de tijd 't') gegeven zijn. % % Link: % http://www.cs.kuleuven.ac.be/~ade/WWW/NW/NW/slinger/slinger/index.html % % Koen Verheyden, 9 augustus 2004 theta = theta_en_omega(1); omega = theta_en_omega(2); g l m c

= = = =

9.81; 1; 1; .1;

% % % %

gravitatieconstante lengte van de slinger massa van de bol dempingsconstante

a = g/l; b = c/(m*l); % rechterlid is een kolomvector met 2 elementen rechterlid = zeros(2,1); % VUL AAN : rechterlid(1,1) = omega ; rechterlid(2,1) = -a*sin(theta)-b*omega;


39



Oplossen van de differentiaalvergelijking op [0, 50] met door middel van de voorwaartse Euler met staplengte

.

>> [t, theta_en_omega] = vweuler('slinger', [0 50], 0.02, [pi/2 5.12]); >> help vweuler VWEULER = voorwaartse euler Gebruik: [t,y] = vweuler('fun',interval,h,y0) Invoer: fun = een m-file die de differentiaalvergelijking y'(t) = fun(t,y(t)) voorstelt de uitvoer is een kolomvector!!! interval = het interval waarover je wil integreren: [a b] h = de staplengte y0 = de beginwaarden (als vector!!!) Uitvoer: t = de punten waarin de oplossing werd berekend. y = de oplossing als matrix (elke kolom bevat de oplossing van een variabele)

De plot van het resultaat in het fasevlak ( op de horizontale as en as).

op de verticale

>> figure; plot(theta_en_omega(:,1), theta_en_omega(:,2), 'g')

Schommelingen, geen convergentie naar nul.


40



Oplossen van de differentiaalvergelijking met andere methoden: Heun (RungeKutta(2) en Runge-Kutta(4). [trk3, yrk3] = rk3('slinger', [0 50], 0.02, [pi/2 5.12]); [trk2, yrk2] = heun('slinger', [0 50], 0.02, [pi/2 5.12]); [trk1, yrk1] = vweuler('slinger', [0 50], 0.02, [pi/2 5.12]); [trk4, yrk4] = rk4('slinger', [0 50], 0.02, [pi/2 5.12]); De functie visslinger geeft de beweging van de slinger weer in functie van de tijd

>> visslinger(0.02, yrk4(:,1)) >> visslinger(0.02, yrk2(:,1)) >> visslinger(0.02, yrk1(:,1)) Runge-Kutta(4) geeft de meest realistische slingerbeweging weer. Bij Heun zijn er twee doorslagen door het hoogste punt, wat fysisch nog mogelijk kan zijn. Voorwaartse Euler daarentegen geeft meerdere doorslagen die fysisch onmogelijk zijn, de demping komt niet duidelijk naar boven. 

Welke functies zijn er in Matlab om differentiaalvergelijkingn op te lossen? Welke raadt Matlab aan? Differential equation solvers. Initial value problem solvers for ODEs. (If unsure about stiffness, try ODE45 first, then ODE15S.) ode45 - Solve non-stiff differential equations, medium order method. ode23 - Solve non-stiff differential equations, low order method. ode113 - Solve non-stiff differential equations, variable order method. ode23t - Solve moderately stiff ODEs and DAEs Index 1, trapezoidal rule. ode15s - Solve stiff ODEs and DAEs Index 1, variable order method. ode23s - Solve stiff differential equations, low order method. ode23tb - Solve stiff differential equations, low order method.


41

Oefening 2: Orde van een methode bepalen De orde van een methode kan men bepalen door te kijken door wat de lokale discretisatiefout en de globale discretisatiefout worden begrensd. Wanneer de globale discretisatiefout begrensd wordt door en de lokale discretisatiefout begrensd wordt door , dan is de methode van orde p. Aangezien de figuur op logaritmische schaal g emaakt is, kan men eenvoudig inzien dat:

Dit heeft tot gevolg dat men naarmate men meer stappen neemt, de fout kleiner wordt. Op basis van de figuur en bovenstaande redenering kan ingezien worden dat de methode van voorwaartse Euler van orde 1 is, de methode van Cauchy van orde 2 en de methode van Heun van orde 2 is.


42

Oefening 3: Stabiliteitsgebieden Stabiliteitsgebieden in het kort (zie ook de uitleg in de Maple-werkbladen): Beschouw de eenvoudige testvergelijking: en een oplossingsmethode met staplengte h. Het stabiliteitsgebied zijn alle waarden (in het complexe vlak) waarvoor de berekende oplossing (voor een bepaalde methode) uitsterft (naar nul gaat). In de cursus zijn het stabiliteitsgebied en het stabiliteitsinterval enkel besproken voor eenstapsmethodes (zoals bvb. Runge Kutta methodes). De doorsnede van B met de reële rechte is het stabiliteitsinterval. B wordt bepaald door een benadering voor de testvergelijking op te stellen zodat yk+1=ϕ(hλ)yk. Runge Kutta: eerste orde en tweede orde: [-2,0] Runge Kutta: derde orde: [-2.5, 0] Achterwaartse Euler: [0, 2] Trapeziumregel: enkel voor 0: geen convergentie, ook niet in 0 aan de rand.

Stabiliteit voor Runge Kutta 1,2,3 en 4


43

Oefenzitting 6: Het iteratief oplossen van een niet-lineaire vergelijking Oefening 1: Numerieke experimenten met bissectie, secantmethode, Newton-Raphson-methode -Algoritmen en hun implementatie in Matlab: function [x,xList]=bisec(f,a,b,eps) k = 0; xList = []; fa = feval(f,a); fb = feval(f,b); if( sign(fa) == sign(fb) ) error('f(a)*f(b) is not negative') end while abs(b-a) > 2*eps % bisection update x = (a + b)/2; fx = feval(f,x); if sign(fx) == sign(fb) b = x; fb = fx; else a = x; fa = fx; end if (mod(k,10) == 0) I = b-a end % loop postamble k = k + 1; xList(k) = x; end x = (a+b)/2; xList(k+1) = x;


44

function [b,bList]=secant(f,a,b,eps, maxIt) k=2; bList = [a,b]; while abs(b-a) > eps % secant update c = a; a = b; b = b + (b - c)/(feval(f,c)/feval(f,b)1); % loop post-amble k = k + 1; bList(k) = b; if(k>maxIt) warning('max iterations reached'); break; end end function [x,xList]=newton(f,df,a,eps,maxIt) k = 1; x=a; xList(1) = x; while true %newton update xprev = x; x = x - feval(f,x)/feval(df,x);

% loop post-amble k = k + 1; xList(k) = x; if( abs(x - xprev) <= eps) break; end if(k>maxIt) warning('max iterations reached'); break; end end


45

-Berekenen van een nulpunt van de functie

:

>> f = @(x)exp(x)-1.5; >> log(1.5) ans = 0.4055

 Startwaarden bepalen: a links van de wortel en b rechts van de wortel. Bijgevolg is volgende keuze van startwaarden een goede keuze:  Eps is de nauwkeurigheid en kan beschreven worden door een negatieve macht van 10. Stel dat eps = Dan:  maxIt geeft het maximaal aantal iteraties weer. Stel bv 100. Dit levert onderstaande code op voor de bissectiemethode >> [b1, b2]=bisec(f,0.3,0.5,10^-16) b1 = 0.4055 b2 = Columns 1 through 9 0.4000 0.4051

0.4500

0.4250

0.4125

0.4063

0.4031

0.4047

0.4055

0.4054

0.4054

0.4055

0.4055

0.4055

0.4055

0.4055

0.4055

0.4055

0.4055

0.4055

0.4055

0.4055

0.4055

0.4055

0.4055

0.4055

0.4055

0.4055

0.4055

0.4055

0.4055

0.4055

0.4055

0.4055

0.4055

0.4055

0.4055

Columns 10 through 18 0.4053 0.4055

0.4054


0.4055


0.4055


0.4055

Columns 46 through 51 0.4055

0.4055


46

Secant: [s1 s2]=secant(f,0.3,0.5,10^-16,10^2) s1 = 0.4055 s2 = 0.3000 0.4055

0.5000

0.4005

0.4052

0.4055

0.4055

0.4055

Newton-Rhapson: [n1 n2]=newton(f,@(x) exp(x),0.3,10^-16,10^2) n1 = 0.4055 n2 = 0.3000

0.4112

0.4055

0.4055

0.4055

0.4055

Er is steeds sprake van convergentie naar 0.4055. De bissectiemethode en Secantmethode vertonen lineaire convergentie, Newton-Rhapson daarentegen vertoont kwadratische convergentie. Dit kan men eenvoudig zien door de fout te plotten. Dit doet men als volgt: Newton-Rhapson

Secant

Bissectie

format long (n2-n1)'

format long (s2-s1)'

format long (b2-b1)'

ans =

ans =

ans =

-0.105465108108164 0.005762222914412 0.000016569764955 0.000000000137278 0 0

-0.105465108108164 0.094534891891836 -0.004989999354450 -0.000232335841153 0.000000580182450 -0.000000000067401 0 0

-0.005465108108164 0.044534891891836 0.019534891891836 0.007034891891836 0.000784891891836 -0.002340108108164 -0.000777608108165 0.000003641891835 -0.000386983108165

-0.000191670608164


47

-Bereken een nulpunt van na wat er gebeurt voor startwaarde 1.

door middel van Newton-Rhapson en ga

Startwaarde 1: divergentie: sprongen tussen 1 en 2 [n1 n2]=newton(g,@(x) 3*x^2-6*x+1,1,10^-16,10^2) n1 = 1 n2 = Columns 1 through 16 1 2

2 1

1 2

2

1

2

1

2

1

2

1

2

1

1

2

1

2

1

2

1

2

1

1

2

1

2

1

2

1

2

1

1

2

1

2

1

2

1

2

1

1

2

1

2

1

2

1

2

1

1

2

1

2

1

2

1

2

1

Columns 17 through 32 1 2

2 1

1 2

2


2 1

1 2

2


2 1

1 2

2


2 1

1 2

2


2 1

1 2

2

Columns 97 through 101 1

2

1

2

1


48

Oefening 2: Consistentie en convergentie van stationaire een-stapssubstitutiemethoden Alle nulpunten van f(x) dienen vaste punten van F(x) te zijn. Bijgevolg dient er bewezen te worden dat als f(x)=0, dan F(x)=x. Dit kan eenvoudig nagegaan worden a.d.h.v. de bijgevoegde maple worksheet, genaamd: substitutiemethoden.mws. Eerst wordt het nulpunt van de functie gezocht m.b.v. Maple: enerzijds grafisch, anderzijds exact.

We laten dit eerst exact uitrekenen door Maple: > solve(f(x)=0, x); De oplossing van deze vergelijking kan uitgedrukt worden door middel van een speciale functie LambertW (geinteresseerd? zie http://www.mathworld.com/LambertsWfunction.html), maar we willen natuurlijk een numeriek resultaat voor deze uitkomst: > xstar := evalf(solve(f(x)=0, x)); > evalf(sqrt(2));


49

-

a)

(

)

Deze substitutiemethode is een voorbeeld van directe substitutie, aangezien x uit de originele vergelijking kan opgelost worden a.d.h.v. een vorige benadering. Grafisch worden de vaste punten gezocht. Een punt is een vast punt wanneer geldt dat:

We zoeken de vaste punten (we willen dat deze overeenkomen met de nulpunten van f(x)): > solve(F(x)=x, x); Volledigheidshalve bekijken we het punt 0 nog eens apart, Maple heeft dat hierboven buiten beschouwing gelaten omdat we anders log(0) nemen, daarom doen we dit via de limiet (als de uitkomst 0 is dan wordt 0 afgebeeld op 0 en is 0 dus ook een vast punt): > limit(F(x), x=0);

De enige vaste punten zijn {LambertW(1)}, dus F(x) is volledig consistent. We gaan nu de convergentie na (zie pagina 32 en verder): > D(F)(xstar); We verwachten dus dat we divergeren, we gaan deze substitutieformule nu gebruiken om 10 iteratiepunten te berekenen. Als startwaarde nemen we x[0] = 0.9: > x[0] := 0.9; for i from 1 to 10 do x[i] := evalf( F( x[i-1] ) ); od;


50

>

> p2 := plot([[x[0], 0], [x[0], x[1]]], color=red): p3 := plot([seq( [[x[k-1], x[k]], [x[k],x[k]], [x[k],x[k+1]] ], k=1..9 )], color=red): display(p1, p2, p3);


51

-

d)

√

We zien op de grafiek dat we 1 vast punt hebben. We gaan nu onderzoeken of dit vast punt inderdaad de oplossing is van onze vergelijking en of we via deze substitutieformule naar deze oplossing zullen convergeren. We zoeken de vaste punten (we willen dat deze overeenkomen met de nulpunten van f(x)): > solve(F(x)=x, x); Volledigheidshalve bekijken we het punt 0 nog eens apart, Maple heeft dat hierboven buiten beschouwing gelaten omdat we anders log(0) nemen, daarom doen we dit via de limiet (als de uitkomst 0 is dan wordt 0 afgebeeld op 0 en is 0 dus ook een vast punt): > limit(F(x), x=0);

De enige vaste punten zijn {0,LambertW(1)}, dus F(x) is volledig consistent. We gaan nu de convergentie na (zie pagina 32 en verder): > D(F)(xstar); We verwachten dus dat we convergeren, we gaan deze substitutieformule nu gebruiken om 10 iteratiepunten te berekenen. Als startwaarde nemen we x[0] = 0.9: > x[0] := 0.9; for i from 1 to 10 do x[i] := evalf( F( x[i-1] ) ); od; >


52

> p2 := plot([[x[0], 0], [x[0], x[1]]], color=red): p3 := plot([seq( [[x[k-1], x[k]], [x[k],x[k]], [x[k],x[k+1]] ], k=1..9 )], color=red): display(p1, p2, p3);

Conclusie convergent en consistent (niet volledig consistent, want niet reciproque consistent).


53

Oefening 3: Waar of niet waar

1. Waar: zie p65 formule 2.1: (

)

2. Niet waar: Het aantal juiste cijfers wordt met twee vermenigvuldigd. 3. Niet waar: Dit geldt niet vanaf de eerste stap, aangezien er schommelingen zijn. Dit geldt echter wel wanneer men steeds dichter bij de oplossing komt.


54

Oefening 4: Newton-Raphson methodes voor het berekenen van √ -

a)

-

> solve(f(x)=0, x);

> unassign('x'); unassign('a'); > f := x -> x^2 - a; -

-

We stellen nu de substitutieformule op voor Newton-Raphson: > F := unapply(simplify(x - f(x)/D(f)(x)), x); -

-

De vaste punten: > solve(F(x)=x, x); -

-

> solve(F(F(x))=x, x); -

-

Voor 2 > a := 2; -

-

> x[0] := sqrt(3*a); for i from 1 to 10 do x[i] := evalf( F( x[i-1] ) ); od; -


55

-

> p1 := plot([x, F(x)], x=-3..3, y=-3..3, color=[black, blue]): p2 := plot([[x[0], 0], [x[0], x[1]]], color=red): p3 := plot([seq( [[x[k-1], x[k]], [x[k],x[k]], [x[k],x[k+1]] ], k=1..9 )], color=red): display(p1, p2, p3);


56

-

b)

> unassign('x'); unassign('a'); > f := x -> 1 - a/x^2;

> solve(f(x)=0, x); We stellen nu de substitutieformule op voor Newton-Raphson: > F := unapply(simplify(x - f(x)/D(f)(x)), x);

De vaste punten: > solve(F(x)=x, x); > solve(F(F(x))=x, x);

Voor 2 > a := 10; > x[0] := sqrt(3*a); for i from 1 to 10 do x[i] := evalf( F( x[i-1] ) ); od;


57



58

-

c)

> unassign('x'); unassign('a'); > f := x -> x- a/x;

> solve(f(x)=0, x); We stellen nu de substitutieformule op voor Newton-Raphson: > F := unapply(simplify(x - f(x)/D(f)(x)), x);

De vaste punten: > solve(F(x)=x, x); > solve(F(F(x))=x, x); Voor 2 > a := 10; > x[0] := sqrt(3*a); for i from 1 to 24 do x[i] := evalf( F( x[i-1] ) ); od;


59



60

Extra 1: Regula falsi function [nulpunt] = regula(f,a,b,eps,Kmax) %REGULA % Deze functie beschrijft het algoritme voor de regula falsi %Input - f is de functie die wordt ingegeven als een String 'f' % - a en b zijn de linker- en rechterstartwaarde % - eps is de tolerantie voor de functie in het nulpunt % - Kmax is het maximum aantal iteraties %Output - nulpunt is het nulpunt, bekomen door de methode van de regula falsi. k=1; fa = feval(f,a); fb = feval(f,b); if( sign(fa) == sign(fb) ) error('f(a)*f(b) is not negative') end for k=1:Kmax %bissectiestap indien het verschil in startwaarden kleiner is dan de %nauwkeurigheid die vooropgesteld werd.Het algoritme wordt meteen %beëindigd. if abs(b-a)<2*eps nulpunt = (b+a)/2; k=Kmax; %Indien dit niet het geval is wordt een secantstap toegepast. else nulpunt = b-(b-a)*fb/(fb-fa); fnulpunt = feval(f,nulpunt); %Secantvoorwaarde if fa*fnulpunt<0 b=nulpunt; fb=fnulpunt; else a=nulpunt; fa=fnulpunt; end end end end

Toegepast op Oefening 1: >> regula(f,0.3,0.5,10^-10,11) ans = 0.405465108108164


61

Oefenzitting 7: Het iteratief oplossen van stelsels (niet-) lineaire vergelijkingen Oefening 1: Conditie van een wortel van een niet-lineaire vergelijking a)Stel dat f(x)=0, twee enkelvoudige wortels, x* en x(2) heeft die dicht bij elkaar liggen. D.w.z. dat klein. Uit formule (2.27) kan je afleiden dat normaal gezien de wortels slecht geconditioneerd zijn, d.w.z. dat groot zal zijn t.o.v. .

Aangezien de afgeleide van f(x) klein is wanneer twee enkelvoudige wortels dicht bij elkaar liggen, kan men in het limietgeval het volgende beschouwen:

Uit formule (2.27) volgt dan dat . Dit heeft tot gevolg dat waaruit men kan concluderen dat de wortels slecht geconditioneerd zijn.

,

b) Gebruik (2.27) om volgende resultaten te verklaren:

Het verband tussen

en ̅

kan geschreven worden als: ̅

Dit ingevuld in (2.27) geeft resultaten in de grootteorde van het werkelijke verschil:

-0.0417

0.6667

-2.2500


2.6667

-1.0417

62

Oefening 2: Stelsels niet-lineaire vergelijkingen Beschouw het stelsel Leid formules af voor volgende iteratieve methodes:

a) elementaire substitutiemethode (vertaal naar vast-punt iteratie) Maak hierin een onderscheid tussen totale en enkelvoudige stapmethodes. (je mag hier uiteraard ook de volgorde van de vergelijkingen veranderen!) De totale stapmethode brengt eerst alle correcties aan op de componenten van x en brengt deze gelijktijdig in rekening. Dit levert onderstaande iteratieformules op:

De enkelvoudige stapmethode berekent eerst de x-waarde in functie van de waarden bekomen in de vorige iteratiestap. Vervolgens wordt de y-waarde bepaald door de zonet bepaalde x-waarde.

b) vereenvoudigde Newton-Raphson Maak weer het onderscheid tussen totale en enkelvoudige stapmethodes. Op p. 125 van de cursus staat de methode uitgewerkt voor de vereenvoudigde Newton-Raphson methode. Eerst worden de twee vergelijkingen geschreven als: { Vervolgens worden de formules (3.14a) en (3.14b) uitgewerkt voor de totale stap methode. Dit levert onderstaande iteratieformules op:


63

Voor de enkelvoudige stap methode worden de formules (3.15a) en (3.15b) uitgewerkt:

c) Newton-Raphson De methode van Newton Raphson steunt op het oplossen van onderstaand stelsel, met als coëfficiëntenmatrix de Jacobiaan van f.

Dit levert als oplossingen voor dit stelsel:

De iteratieformules worden dan bepaald door onderstaande uitdrukking in te vullen:


64

Oefening 3: Convergentie voor de methode van Gauss-Seidel voor stelsels lineaire vergelijkingen Voor het oplossen van stelsels lineaire vergelijkingen kunnen iteratieve methoden zoals Gauss-Seidel efficiënt zijn indien de matrix ijl (spaars) is, zoals bij stelsels afkomstig van discretisaties van partiële differentiaalvergelijkingen. In deze oefening bekijken we de rekenkost van de methode van Gauss-Seidel voor het oplossen van het stelsel afkomstig van de discretisatie van de tweedimensionale Laplacevergelijking. De rekenkost is gelijk aan de rekenkost per iteratie maal het aantal iteraties om een bepaalde nauwkeurigheid te bereiken, stijgt. (http://people.cs.kuleuven.be/~adhemar.bultheel/WWW/NW/NW/laplace/laplac/index.html) De convergentiefactor van de methode van Gauss-Seidel voor het beschouwde stelsel is

a) Hoeveel iteratiestappen zijn er nodig om de fout met een factor 104 te verkleinen? Op basis van de definitie van de convergentiefactor, kan de iteratiefout geschreven worden in functie van de convergentiefactor en de iteratiefout bij de nulde iteratiestap.

Dit levert op Gauss-Seidel toegepast onderstaande uitdrukking op voor het aantal iteratiestappen:

(

)

b) Wat is de convergentiefactor en hoeveel stappen zijn er nodig om de fout met een factor 10-4 te verkleinen indien n=100?

c) Wat is de rekenkost van 1 iteratiestap in functie van het aantal discretisatiepunten n in de x- en y-richting? Aangezien men te maken heeft met een pentadiagonale matrix, met n-1 rijen en n-1 kolommen waarover geïtereerd dient te worden is de rekentijd van O(n²). d) Stel dat de rekentijd van de Gauss-Seidelmethode voor n=100 gelijk is aan k1 sec, wat is dan de rekentijd voor n=200?


65

Oefenzitting 8: Eigenwaarden, Optimalisatie en Herhaling Oefening 1: Het berekenen van eigenwaarden Waarom werkt de methode van de machten niet voor onderstaand geval? * Bereken

+

* +

en interpreteer.

Men kan eenvoudigweg zien dat de eigenwaarden 3 en 2 zijn, aangezien men te maken heeft met een bovendriehoeksmatrix, waarbij de eigenwaarden gelijk zijn aan de diagonaalelementen. Wanneer men de eigenwaarden manueel uitwerkt bekomt men volgende eigenwaarden: |

|

De eigenwaarden die men zou moeten bekomen zijn dus 3 en 2. Vervolgens worden de eigenvectoren bepaald. Deze bekomt men door de stelsels die men bekomt door de eigenwaarden in te vullen op te lossen:

Bijgevolg kan de vector

*

+

* +

* +

*

+

* +

* +

geschreven worden als:

De methode van de machten of de methode van von Mises berekent een dominante eigenwaarde . Dit betekent dat de grootste eingenwaarde niet complex mag zijn. Dit betekent echter niet dat alle eigenwaarden verschillen. Een willekeurige vector is te schrijven als een lineaire combinatie van de eigenvectoren vermits deze een vectorruimte opspannen. De methode van de machten geldt enkel als de coëfficiënt bij de eigenvector van de grootste eigenwaarde verschillend is van 0. Dit is echter niet het geval, zodat de methode van de machten niet werkt voor onderstaand geval. Indien men deze eenvoudige berekening niet zou gemaakt hebben en klakkeloos de methode van de machten zou toepassen, bekomt men onderstaande uitdrukking voor . ∑ Hieruit zou men dan impliceren dat wat men op het zicht ziet.

* +

* +

* +

de grootste eigenwaarde is, wat niet strookt met


66

Oefening 2: Optimalisatie: waar of niet waar? -

Bij het minimaliseren van een unimodale functie in één veranderlijke met de methode van de gulden snede wordt altijd het punt met de hoogste functiewaarde weggelaten. Niet waar: Er moet steeds getest worden dat f(u)>f(v). Er bestaan gevallen waarbij a bv niet de hoogste functiewaarden heeft, maar toch weggelaten wordt. Stel dat men volgende functie heeft: 𝑓(𝑢 𝑘 ) < (𝑣 𝑘 ) Nu kiest men een nieuw interval als [𝑢 𝑘 𝑣 𝑘 ] [𝑎 𝑘 𝑏 𝑘 ] In dit geval wordt de grootste functiewaarde niet weggelaten

-

Zal de methode van de gulden snede convergeren voor een monotoon (dalende of stijgende) functie? Indien ja, naar welk punt. Indien niet, waarom? Waar: Er is convergentie naar de eindpunten. Voor een monotoon dalende functie is er convergentie naar b en voor een monotoon stijgende functie is er convergentie naar a.


67

Oefening 3: Herhaling: Stelsels lineaire vergelijkingen Gegeven het stelsel:

waarbij A een matrix is met conditiegetal . We hebben dit probleem opgelost met Matlab en kwamen met de volgende bevindingen. De 2-norm van de residu is gelijk aan en ‖ ‖ . De absolute fout op de oplossing is en ‖ ‖ . Voor de stabiliteit wordt achterwaartse stabiliteit nagegaan. Hierbij wordt gecontroleerd of de verhouding van de norm het residu tegenover de norm van ongeveer gelijk is aan de machinenauwkeurigheid. ‖ ‖ Er is geen voorwaartse stabiliteit aangezien de 2-norm van het residu niet van de grootteorde van het machinegetal is. De methode is stabiel aangezien de relatieve fout en absolute fout op de oplossing zeer kleine zijn. Aangezien het conditiegetal hoog is , betekent dit dat het probleem slecht geconditioneerd is. Conditie is een kenmerk van het probleem, in dit geval het stelsel lineaire vergelijkingen. Men verwacht een verlies van maximaal zes cijfers. Wanneer men dit vergelijkt met de relatieve fout, stelt men een verlies van drie cijfers vast.


68

Oefening 4: Herhaling: Stelsels lineaire vergelijkingen: convergentie van iteratieve methoden We lossen een stelsel van 5 vergelijkingen met 5 onbekenden op met de methode van Jacobi. We stellen de onbekenden voor als een vector Als startwaarde gebruiken we de vector ( ) . De resultaten van de eerste vijf iteraties zijn gegeven in onderstaande tabel.

De vector convergeert naar (25.000000, 35.714285, 42.857143, 35.714285, 25.000000). Wat is de convergentiefactor voor de opeenvolgende benaderingen voor ‖ ‖ ? Hoeveel iteraties denk je nodig te hebben om de absolute fout op kleiner te maken dan ? De methode van Jacobi is een totale stapmethode, waarbij de j-de vergelijking van volgende vorm is:

Bijgevolg kan de iteratie herschreven worden als:

In bovenstaande uitdrukking is het residu gelijk aan: (

)

De methode van Jacobi kan ook beschreven worden in matrixnotatie. Dit levert onderstaande uitdrukking op:

De k-de iteratiefout voor de methode van Jacobi: (

)

De definitie van de convergentiefactor geeft onderstaande uitdrukking:


69

Toegepast op

levert dit het volgende:

Hieruit kan men concluderen dat men per twee iteraties de iteratiefout verandert als:

Het aantal iteraties nodig om de absolute fout op worden door onderstaande uitdrukking:

kleiner te maken dan

, kan bepaald

<

Hieruit volgt dat er 13 iteraties nodig zijn.


70

Oefening 5: Herhaling: Stabiliteit van de methode van Euler Stel we hebben een stelsel differentiaalvergelijkingen:

met beginwaarde , waarbij y een n-vector is en A een reëel symmetrische matrix. Het stabiliteitsgebied van de methode van Euler bepaalt hoe groot de stapgrootte mag gekozen worden. De volgende figuren tonen ‖ ‖ in functie van voor enkele stapgrootten . De eigenwaarden van de matrix zijn negatief, dus moet de oplossing naar nul convergeren.

Op basis van deze figuren kan je besluiten in welk interval de eigenwaarden van A moeten liggen. Welk interval is dat? Motiveer. Aangezien men te maken heeft met een reële symmetrische eigenwaarden. Het stabiliteitsinterval voor de methode van Euler is: ] Bijgevolg kan geëist worden dat

, heeft men reële

[

in dit interval ligt: <

<

Op basis van bovenstaande figuren wordt nagegaan welke stapgrootte gekozen dient te worden opdat de vergelijking uitdooft. Er wordt begonnen met de hoogste stapgrootte. Zo bekomt men h=0.018. Dit levert onderstaand interval op waarbinnen de eigenwaarden van de matrix A moeten liggen: <

< <

<


71

Examenvraag 1: LU-factorisatie Stel dat we volle stelsels van n lineaire vergelijkingen met n onbekenden willen oplossen met LU-factorisatie of Gausseliminatie. 1. Veronderstel dat een PC 0.14 seconden nodig heeft voor de LU-factoristie voor n=800. Hoeveel tijd zal die PC nodig hebben om stelsels met n=1600 en n=8000 onbekenden op te lossen? Zeg ook waarom. Het oplossen van een stelsel met n onbekenden en n vergelijkingen heeft een tijdsduur die volledig afhangt van die van de LU-factorisatie. De LU-factorisatie heeft een rekenkost van O(n³). Dit kan toegepast worden op verschillende waarden van n. Aangezien men dus weet dat voor n=800, t=0.14 seconden, dan kan men 1600 schrijven als 2n en 8000 schrijven als 10n. Dit levert onderstaande uitdrukkingen op:

2. Indien het bijhouden van de matrix in de PC ongeveer 5.12 MB vergt, hoeveel MB zijn er dan nodig voor n=1600 en n=8000? Zeg ook waarom. In een - matrix zijn elementen aanwezig. Aangezien voor een matrix het benodigde geheugen ongeveer 5.12 MB bevat, kan er geschreven worden dat:

Voor n=1600 levert dit het volgende op:

3. Indien de PC beschikt over 2000 MB geheugen, wat is de maximale grootte van een matrix die volledig in het geheugen kan bewaard worden? Zeg ook waarom. Het benodigd geheugen per element wordt gegeven door onderstaande uitdrukking en is gelijk aan:

Het geheugen kan bepaald worden door:


72

Aangezien het geheugen gekend is en het aantal elementen gezocht wordt, kan dit geschreven worden als volgt:

Stel dat de matrix een elementen gelijk aan:

-matrix is met aantal elementen

. Dan is het aantal

Indien men dan een vierkante matrix wil opslaan in de PC, kan het aantal rijen en kolommen van de matrix bepaald worden door:

√ Bijgevolg kan er maximaal een PC.


opgeslagen worden op de

73

Examenvraag 2: Juiste cijfers Onderstaande tabel geeft het resultaat van de volgende berekeningen in Matlab: for i=1:8 x(i)=10^(-i) w(i)=(1-cos(x(i)))/(x(i)^2) end i 1 2 3 4 5 6 7 8

x(i) 1.00e 1.00e 1.00e 1.00e 1.00e 1.00e 1.00e 1.00e

− − − − − − − −

1 2 3 4 5 6 7 8

w(i) 4.995834721974289e 4.999958333473664e 4.999999583255033e 4.999999969612645e 5.000000413701856e 5.000444502911705e 4.996003610813205e 0

− − − − − − −

01 01 01 01 01 01 01

Vermits verwachten we dat de waarden voor stijgende dichter bij 0.5 liggen. Als je als benadering voor 0.5 beschouwt: - Hoeveel juiste cijfers heeft ? - Hoeveel juist cijfers heeft ? - Verklaar waarom de meest nauwkeurige benadering geeft. Het aantal juiste cijfers wordt bepaald door de definitie van een juist cijfer toe te passen. Als ̅ het laatste juiste cijfer is (dus ̅ het eerste foutieve), dan zal gijgevolg gelden dat:

Als ̅ dus voldoet aan: < Dan zijn er ̅

juiste cijfers.

Toegepast op dit geval geldt er: < Voor

<

̅

zijn er 8 juiste cijfers. <

<

̅

Voor w(5) zijn er 7 juiste cijfers.


74

De verklaring hiervoor wordt onderstaande Matlab-code en de foutencurve die ermee werd getekend: function [] = examenvraag2( input_args ) %examenvraag2 Summary of this function goes here % Detailed explanation goes here % Deze methode tekent de grafiek van de absolute fout op een semilogplot. % % lim (x--> 0) (1-cos(x))/x^2 = 0.5 % Bovenstaande uitdrukking wordt getest met opeenvolgende machten (1ste % macht t.e.m. 20ste macht) van 10 for i = 1:20 x = 10.0^-i; absolutefout(i) = (1-cos(x))/x^2; end; % plotten van de grafiek van de absolute fout figure(1) semilogy(abs(absolutefout-0.5)) title('absolute fout') end

Bespreking foutencurve: In het begin van de grafiek is er enkel de benaderingsfout die een rol speelt. Deze neemt lineair af tot men komt bij . Vanaf zijn de afrondingfouten groter dan de benaderingsfouten. Bij of groter telt Matlab een getal op bij 1 dat kleiner is dan de machinenauwkeurigheid (

). Bijgevolg zal Matlab dit getal afronden naar 1, zodat de

fout constant blijft.


75

Examenvraag 3: waar of niet waar? 1. De integratieregel van de methode van Simpson integreert een veelterm van graad twee exact. Waar of niet waar? Geef een korte verklaring. De regel van Simpson is een samengestelde regel waarbij in elk deelinterval een interpolerende veelterm van de tweede graad gebruikt wordt. Daarvoor heeft men in elk deelinterval drie interpolatiepunten nodig. De uitspraak is niet waar in de zin, dat de methode van Simpson enkel een veelterm van graad twee exact zal integreren wanneer de vierde afgeleide bestaat, aangezien de intebgratiefout naar nul gaat. Details omtrent deze berekening zijn te vinden in de cursus op p. 42 en 43. 2. Het bepalen van een meervoudige wortel van een niet-lineaire vergelijking is slecht geconditioneerd. Waar of niet waar? Geef een korte verklaring. Deze stelling is waar. Het komt erop neer dat wanneer de absolute waarde van de afgeleide klein is of nul, het probleem slecht geconditioneerd is. Daarentegen wanneer de absolute waarde van de afgeleide groot is, zal het probleem goed geconditioneerd zijn. Aangezien men te maken heeft met het bepalen van een meervoudige wortel, kan men concluderen dat de afgeleide steeds nul is, wat betekent dat het probleem slecht geconditioneerd is. Voor een meer gedetailleerde uitleg wordt er verwezen naar de cursus op p 101 en 102.


76

Examenvraag 4: Stabiliteit van methoden voor het oplossen van differentiaalvergelijkingen 1. Beschouw het stelsel lineaire differentiaalvergelijkingen

*

+

Wat is de maximale stapgrootte zodat de volgende methoden stabiel zijn? (a) Runge-Kutta(4)-methode (b) Methode van de trapeziumregel Oplossing: De eigenwaarden van de matrix A zijn -100 en -2. Deze kan men snel bepalen, aangezien de matrix A een bovendriehoeksmatrix is, waarvoor geldt dat de eigenwaarden gelijk zijn aan de diagonaalelementen. De oplossing van dit stelsel differentiaalvergelijkingen kan eenvoudig bepaald worden en is de volgende:

De maximale stapgrootte kan bepaald worden door voor elk van beide methodes het stabiliteitsinterval op te stellen en vervolgens na te gaan in welk interval en voor welke eigenwaarde h maximaal is. (a) Voor Runge-Kutta(4) kan het onderstaande geschreven worden: Het stabiliteitsinterval voor Runge-Kutta(4) (Tabel 7.4) is gelijk aan [-2.78, 0[, Bijgevolg is de maximale stapgrootte gelijk aan 1.39 opdat de methode stabiel is. (b) Voor de methode van de trapeziumregel is het stabiliteitsinterval gelijk aan: het linkerhalfvlak van het complexe vlak: ] [ De maximale stapgrootte is in dat geval gelijk aan oneindig opdat de methode stabiel is.


77

2. Beschouw het stelsel lineaire differentiaalvergelijkingen

*

+

De eigenwaarden van A zijn . Wat is de maximale stapgrootte opdat de voorwaartse methode van Euler stabiel is? Het stabiliteitsinterval voor voorwaartse Euler is gelijk aan : ]

[

Bijgevolg is de maximale stapgrootte gelijk aan 16.


78

Examenvraag 5: Stelsels differentiaalvergelijkingen (Nog te maken) Beschouw het volgend stelsel differentiaalvergelijkingen:

Waarbij C de symmetrische tridiagonale -matrix is met op de hoofddiagonaal een op de eerste diagonalen in boven- en onderdriehoek. We lossen deze differentiaalvergelijking op met de achterwaartse methode van Euler met stapgrootte en met startwaarde . In iedere stap van deze impliciete methode moeten we een lineair stelsel oplossen. Indien de stapgrootte verkleint, wordt de methode van Euler nauwkeuriger, maar er moeten meer stelsels opgelost worden. Stel dat we voor het oplossen van de stelsels de iteratieve methode van Jacobi gebruiken. 1. Schrijf de formule van de achterwaartse methode van Euler voor het stelsel differentiaalvergelijkingen: Algemeen kan de formule voor de achterwaartse methode van Euler geschreven worden als:

Toegepast op dit stelsel differentiaalvergelijkingen:

2. Schrijf het stelsel neer dat in iedere stap moet worden opgelost. Schrijf de matrix van het stelsel volledig uit. Uit vraag 1 volgt het stelsel dat men dient op te lossen:

Er dient opgemerkt te worden dat de coëfficiëntenmatrix van dit stelsel als volgt kan geschreven worden:

[


]

79

3. Een voldoende voorwaarde voor convergentie van de methode van Jacobi, is dat de matrix van het op te lossen stelsel diagonaaldominant is. Voor welke waarden van is de matrix van het op te lossen stelsel diagonaaldominant? Om te bepalen of de matrix van het stelsel diagonaaldominant is dient er nagegaan te worden of het diagonaalelement groter is dan de som van de andere elementen in een rij. Dit levert onderstaande ongelijkheid op:

Bijgevolg vindt men onderstaande voorwaarde voor de stapgrootte h, deze moet in elk geval groter zijn dan nul: < 4. De norm van de iteratiematrix van h is:

<

van de methode van Jacobi voor dit stelsel als functie

waarbij het aantal rijen en kolommen is van de matrix. Maak een tabel met vier kolommen. In de eerste kolom plaats je de volgende waarden van : , en , in de tweede kolom van de overeenkomstige (bovengrens voor de) norm van de iteratiematrix. In de derde kolom plaats je het aantal Jacobiiteraties dat je verwacht nodig te hebben opdat de oplossing van één stelsel 5 beduidende cijfers meer heeft dan de startwaarde. In de vierde kolom plaats je het totaal aantal Jacobi-iteraties voor alle stappen van de achterwaartse methode van Euler om het stelsel differentiaalvergelijkingen te itereren van x=0 tot x=10 Interpreteer de tabel: wordt de methode van Jacobi goedkoper of duurder voor een stelsel, wanneer kleiner wordt? Wordt achterwaartse Euler met Jacobi-iteratie goedkoper of duurder bij kleiner wordende ? Stapgrootte h

Bovengrens ‖ ‖

1/8

1/3

# Jacobi-iteraties (5 beduidende cijfers meer dan startwaarde) 10.48 11

1/16

1/7

5.92 6

1/32

1/15

4.25 5


Totaal # Jacobiiteraties voor alle stappen Euler (x=0 tot x=10).

80

Om het aantal Jacobi-iteraties om vijf beduidende cijfers meer dan de startwaarde te hebben, dient men onderstaand verband tussen de norm van de iteratiematrix en de iteratiefout toe te passen: ‖

‖

‖ ‖ ‖

‖

Aangezien er vijf meer beduidende cijfers nodig zijn, kan men de norm van k-de iteratiefout schrijven als 1/100000ste van de oorspronkelijke iteratiefout. Dit levert onderstaande vergelijking op voor het bepalen van het aantal nodige Jacobi-iteraties. ‖ ‖ Het aantal Jacobi-iteraties nodig om alle stappen door te voeren die nodig zijn voor de achterwaartse methode van Euler om het stelsel differentiaalvergelijkingen te itereren van x=0 tot x=10, kan als volgt bepaald worden. Er dienen een bepaald aantal stelsels opgelost te worden d.m.v. de methode van Jacobi. Dit aantal stelsels is gelijk aan k, dus 10 in dit geval. Vervolgens kan nagegaan worden wat het verband is tussen het aantal stelsels dat dient opgelost te worden en het aantal Jacobi-iteraties.


81

Examenvraag 6: Optimalisatie in meerdere veranderlijken Beschouw de symmetrische matrix

. Beschouw de functie in

veranderlijken:

De gradiënt is

Als

een eigenvector is van , is 1. Toon aan dat

de bijbehorende eigenwaarde.

als en slechts als

een eigenvector is.

⇐

Als

een eigenvector is van A, dan geldt:

Beide leden vermenigvuldigen met

Beide leden delen door

:

, levert op dat:

Hieruit volgt uit het gegeven voor

dat:

Wanneer men dan de gradiënt van

berekent, bekomt men dat:

Om de pijl in de andere richting te bewijzen dient men het volgende te doen: Numerieke Wiskunde: Oefeningen Philippe Nimmegeers

82

Door beide leden te vermenigvuldigen met

2. Wat is

En

, geldt het volgende:

?

komt overeen met de kleinste eigenwaarde van A. komt overeen met de grootste eigenwaarde van A. 3. Beschouw de matrix A gemaakt met de volgende MATLAB-code x = orth(randn(100)) ; A = x * diag( - 1./[1:100] ) * x ’ ;

De eigenwaarden van

zijn

(a) Voor de minimalisatie van stap berekend wordt als:

voor

.

beschouwen we de afdalingsmethode waarbij de k-de (

-

)

We hebben twee keuzes voor uitgeprobeerd voor dit probleem: De methode van de steilste afdaling waarbij zodanig gekozen wordt dat in de k-de iteratie ( ( )) minimaal is. De methode waarbij gekozen wordt volgens:

(let op het minteken!) Welke van de twee methodes zou sneller moeten zijn? Motiveer. De volgende tabel geeft voor beide methoden de fout op het minimum voor opeenvolgende iteraties. Met welke methoden (d.w.z. van de gegeven) komen A en B dan overeen? De eerste methode komt overeen met methode B, aangezien de methode van de steilste afdaling optimaal functioneert eens de ideale richting gevonden is. Dit is de richting van de minimale functiewaarde. De tweede methode convergeert in het begin misschien beter, maar er zijn toch meer stappen nodig. Dit komt dus overeen met methode A.


83

(b) Wat zijn de respectievelijke convergentiefactoren voor methoden Aangezien bovenstaande tabel gegeven is, kunnen de convergentiefactoren bepaald worden op basis van deze gegevens:

Zo geldt voor methode A:

Voor de methode B:

(c) Als we in de eerste methode zouden kiezen zoals in methode 2, dan komt deze methode overeen met een methode om eigenwaarden te berekenen. Welke? Verklaar. Deze methode komt overeen met de methode van de machten voor het berekenen van eigenwaarden.


84

Oefenzitting 9: Optimalisatie in meerdere veranderlijken Oefening 1: De methode van de steilste afdaling en methode van Newton Gegeven de volgende functies in 1. 2. Onderzoek de convergentie van de methode van de steilste afdaling en de methode van Newton voor deze functies. Zoek voor elke functie uit of: a) Aan de nodige en voldoende voorwaarden voldaan is voor een minimum b) Wat kan je zeggen over de convergentie van de methode van Newton. c) Wat kan je zeggen van de methode van de steilste afdaling. Verklaar je antwoorden. Oplossing: a) De nodige en voldoende voorwaarden zijn de volgende: Nodige voorwaarde:

Voldoende voorwaarde:

Voor de eerste functie betekent dit het volgende: [

] *

+ *

+

Aangezien bovenstaande matrix een diagonaalmatrix is, zijn de eigenwaarden gelijk aan de diagonaalelementen. Deze matrix is positief definiet en heeft als dubbele eigenwaarde Voor de tweede functie betekent dit: [

] *

+


85

Aangezien de som van de rijen gelijk is aan 400, kan er besloten worden dat 400 een eigenwaarde is van bovenstaande matrix. De tweede eigenwaarden kan bepaald worden door het spoor van de matrix te bepalen: dit is de som van de diagonaalelementen en deze is gelijk aan 404. Bijgevolg is de tweede eigenwaarde 4. b) De convergentie voor de methode van Newton kan als volgt bepaald worden: Indien men geen voorzorgen neemt zal de methode van Newton kunnen convergeren naar om het even welk stationair punt dat niet noodzakelijk een minimum is. De correctievector zal wel wijzen in de richting van afnemende functiewaarden indien ( ) . Zie pp 170 en 171 voor een mogelijke aanpassing aan het algoritme die als voorzorgsmaatregel kan geprogrammeerd worden t.v.v. de correctiestap. Aangezien de Hessiaanmatrix niet singulier is, dus een determinant hebben die verschillend van nul is (waaruit volgt dat deze inverteerbaar is) volgt dat de convergentie kwadratisch is. Aangezien de tweede functie een kwadriek is en de methode van Newton een functie benadert d.m.v. een kwadriek, gebeurt dit in een stap. c) De methode van de steilste helling kiest voor daalt: ‖

de richting waarvoor

het snelste

‖

Deze methode is een afdalingsmethode en in elke stap zal er dus een lagere functiewaarde gevonden kunnen worden zolang de gradiënt verschilt van nul. Zie p 168: voor het eerste geval is de grootste eigenwaarde gelijk aan de kleinste (dubbele eigenwaarde 2), zodat men convergentie na één stap heeft (concentrische cirkels). In het tweede geval is de grootste eigenwaarde 100 keer groter dan de kleinste. Dit veronderstelt trage convergentie, die voorgesteld kan worden door elliptische contourlijnen.


86

Oefening 2: Vergelijking methode van Newton en steilste afdaling Onderstaande figuur geeft de richting aan van de methode van Newton en de methode van de steilste afdaling. Vul de legende van de figuur aan: welke richting hoort bij welke methode? Newton: volle pijl, aangezien deze kwadratische convergentie impliceert. Steilste afdaling: halve pijl: tragere convergentie (ellipsen)


87

Oefening 3: Methode der steilste afdaling Beschouw de symmetrische matrix

. Beschouw de functie in

veranderlijken:

De gradiënt is

1. Toon aan dat

als en slechts als

een eigenvector is.

⇐

Als

een eigenvector is van A, dan geldt:

Beide leden vermenigvuldigen met

Beide leden delen door

:

, levert op dat:

Hieruit volgt uit het gegeven voor

dat:

Wanneer men dan de gradiënt van

berekent, bekomt men dat:

Om de pijl in de andere richting te bewijzen dient men het volgende te doen:


88

Door beide leden te vermenigvuldigen met

2. Wat is

En

, geldt het volgende:

?

komt overeen met de kleinste eigenwaarde van A. komt overeen met de grootste eigenwaarde van A. 3. Beschouw de matrix A gemaakt met de volgende MATLAB-code x = orth(randn(100)) ; A = x * diag( - 1./[1:100] ) * x ’ ;

De eigenwaarden van

zijn

(a) Voor de minimalisatie van stap berekend wordt als:

voor

.

beschouwen we de afdalingsmethode waarbij de k-de (

-

)

We hebben twee keuzes voor uitgeprobeerd voor dit probleem: De methode van de steilste afdaling waarbij zodanig gekozen wordt dat in de k-de iteratie ( ( )) minimaal is. De methode waarbij gekozen wordt volgens:

(let op het minteken!) Welke van de twee methodes zou sneller moeten zijn? Motiveer. De volgende tabel geeft voor beide methoden de fout op het minimum voor opeenvolgende iteraties. Met welke methoden (d.w.z. van de gegeven) komen A en B dan overeen? De eerste methode komt overeen met methode B, aangezien de methode van de steilste afdaling optimaal functioneert eens de ideale richting gevonden is. Dit is de richting van de minimale functiewaarde. De tweede methode convergeert in het begin misschien beter, maar er zijn toch meer stappen nodig. Dit komt dus overeen met methode A.


89

(b) Wat zijn de respectievelijke convergentiefactoren voor methoden Aangezien bovenstaande tabel gegeven is, kunnen de convergentiefactoren bepaald worden op basis van deze gegevens:

Zo geldt voor methode A:

Voor de methode B:

(c) Als we in de eerste methode zouden kiezen zoals in methode 2, dan komt deze methode overeen met een methode om eigenwaarden te berekenen. Welke? Verklaar. Deze methode komt overeen met de methode van de machten voor het berekenen van eigenwaarden.


90

Oefening 4: Optimalisatie van een functie in twee veranderlijken. Gegeven de objectfunctie: ∫

met beperking:

Met parameters. De functies een positief definiete Hessiaan.

zijn zodanig dat 1 globaal minimum heeft met

Beginvoorwaarde: Beschrijf hoe je kan bepalen zodanig dat J minimaal is. Zeg ook welke numerieke methoden je hiervoor zou gebruiken. Dit probleem is een optimalisatieprobleem dat kan opgelost worden d.m.v. de methode van de steilste afdaling. Zo kan de gradiënt berekend worden als: ∫ ∫ Dit komt neer op het oplossen van twee integralen. De methode van Simpson kan hiervoor gebruikt worden met een bepaalde stapgrootte h. Eveneens dient er een afgeleide berekend te worden:

Tenslotte is er een stelsel van differentiaalvergelijkingen dat opgelost moet worden d.m.v. de methode van de trapeziumregel of m.b.v. de Runge-kutta(4) methode.

[

]

[


][

]

[

]

91

Numerieke wiskunde: Oefeningen

Recommend Documents