Vanliga frågor

2D1240, Numeriska metoder gk2 för T2, våren 2004

Laboration 1

Fråga: Vad menas med "konvergenshastigheten" i uppgift 2c1?

Svar: En precis förklaringen är följande: Antag att sekvensen x_k uppfyller

lim_(k->oo) |x_(k+1)-a|/|x_k-a|^p = C.

Då säger vi att sekvensen konvergerar med ordning p mot a och att konvergenshastigheten är C.

I fixpunktfallet är p=1 och man får att för stora k, dvs nära roten, gäller

|x_(k+1)-a| ~= C |x_k-a|,

dvs konvergenshastigheten C är den faktor med vilken felet minskar i varje iteration.
Motsvarande för Newtons metod, som konvergerar kvadratiskt med p=2, blir

|x_(k+1)-a| ~= C |x_k-a|^2.

Teoretiska uttryck för C i fixpunktfallet och för Newtons metod finns i Quarteroni/Saleri.


Fråga:
Ska verkligen det experimentella konditionstalet skilja sig så mycket från det teoretiska konditionstalet i uppgift 3?

Svar: Det teoretiska konditionstalet ger den maximala förstärkningen av relativfelet som man kan få om man väljer b och störningen db så illa som möjligt. Det b som är givet i uppgiften är dock inte det värsta möjliga så det experimentella konditionstalet blir därför ganska mycket mindre än det teoretiska. Om man istället väljer b som egenvektorn hörande till matrisens största egenvärde (b~=(0.7926 0.4519 0.3224 0.2522) för n=4) blir situationen mycket värre, speciellt om man sedan väljer db som egenvektorn hörande till matrisens minsta egenvärde (db~=(0.0292 -0.3287 0.7914 -0.5146) för n=4).


Fråga: När jag gör [L,U]=lu(A) blir inte L undertriangulär. Varför inte?

Svar: Det beror på att MATLAB gör LU-faktorisering med radpivotering. Då kommer även en permutationsmatris in i bilden (se kapitel 5.2 i Quarteroni/Saleri) och faktoriseringen som MATLAB räknar ut är PA = LU, där P är permutationsmatrisen. När man anropar lu() som ovan returnerar MATLAB P^T*L istället för L. Skriver man istället

>> [L,U,P]=lu(A);

får man L och U (och P) separat. Några saker att notera:

Fråga: Skall man använda MATLAB på något sätt i uppgift 5a?

Svar: Nej, bara papper och penna.

Laboration 2

Fråga: Vad är den teoretiska stabilitetsgränsen i uppgift 3?

Svar: En numerisk metod för linjära skalära ODEer är absolutstabil för en given steglängd h om de uträknade värdenau_n går mot noll när metoden appliceras på testproblemet

y'=lambda*y,

med lambda < 0. För framåt Euler kan man visa att detta händer precis när

|1+h*lambda| < 1.

Stabilitetsgränsen är det h för vilket detta villkor precis uppfylls. För bakåt Euler kommer u_n gå mot noll för alla h>0 och stabilitetsgränsen är då helt enkelt h=0. (Se QS s 165 och s 167.)

För linjära system av ODEer av typen

y'=A*y,    

där A är en n x n-matris, är framåt Euler absolutstabil om olikheten ovan håller för alla egenvärden till A. Dvs, låt lambda_j, j=1,...n, vara egenvärdena till A. Framåt Euler är då stabil om

|1+h*lambda_j| < 1,    j=1,...,n.

Bakåt Euler är fortfarande stabil för alla h (när alla egenvärdena till A har negativ realdel).

Fråga: Hur verifierar vi stabilitetsgränsen empiriskt?

Svar: Använd framåt Euler för att lösa systemet dels med en steglängd h som ligger lite under och dels med en som ligger lite ovanför stabiltetsgränsen. Plotta de två lösningarna och jämför dem. Den numeriska lösningen med litet hbör hålla sig stabil (begränsad) och den med för storthväxa okontrollerat. För bakåt Euler behöver ni bara visa att lösningen blir begränsad även om ni tar mycket störreh än stabilitetsgränsen för framåt Euler.

Fråga: Varför får jag så konstiga värden på noggrannhetsordningen för framåt Euler på uppgift 4 b?

Svar: Teorin för noggrannhetsordningen är asymptotisk och kräver att felet är litet. Vid tiden t=20 har felet för framåt Euler hunnit växa sig stort även för så små tidssteg som h=2^(-9). Man måste ta mycket mindre tidssteg för att se noggrannhetsordningen, vilket (tyvärr) tar lång tid att köra på datorn. Prova därför istället att jämföra framåt Euler-lösningarna vid tiden t=2. (Problemet uppstår inte med Runge-Kutta eftersom det är en mycket noggrannare metod som ger mindre fel. För Runge-Kutta kan ni därför fortfaranda använda t=20.)

Fråga: Hur väl behöver felen för ode45, Runge-Kutta och framåt Euler "hamna inom samma intervall" i uppgift 4 c?

Svar: Det är svårt att få felen för framåt Euler att ligga i samma intervall som de övriga, eftersom det kräver väldigt små tidssteg och tar lång tid. Försök att få med i alla fall ett värde där felet för framåt Euler är lite mindre än felet för Runge-Kutta med h=0.2. Det kan också eventuellt vara intressant att köra Runge-Kutta med lite mindre tidssteg än h=0.01 för jämförelsen med ode45.


Laboration 3

4.3 Motordrivna inversa pendeln: Plotta pendeln som ett streck och visa hur den rör sig i tiden. (rkpendel.m ingår inte.)

4.5 Strömkretsen: Jämför gärna er Fourieranalys med resultatet från fft, Matlabs kommando för snabb fouriertransform.

4.11 Naturen: Välj tidssteget dt för RK4 så att felet  i T1 är mindre än en dag. (Detta dt kan vara större eller mindre än 1/365.)

4.14 Flödespaketet: Richardsonextrapolation beskrivs i Uppgift 4.10, s 101 i Quarteroni/Saleri.

4.16 Satelliten: Banorna är väldigt störningskänsliga så defaultnoggrannheten i ode45 räcker inte riktigt till för att ge ett bra svar. Sätt upp den tex genom att först sätta "opt=odeset('RelTol',1e-9,'AbsTol',1e-9);" och sedan anropa ode45 med "opt" som extra argument. (Se hjälpen till ode45, odeset.) Dessa komplicerade satellitbanor kallas Arenstorf-banor och upptäcktes så sent som 1962.

4.19 Dubbelpendeln: Plotta dubbelpendeln som två streck och visa hur den rör sig i tiden. (rkpendel.m ingår inte.)