1. No category

Teknisk Beräkningsvetenskap I
Tema 3: Styvhetsmodellering av mjuk mark med
icke-linjära ekvationer
Eddie Wadbro
18 november, 2015
Eddie Wadbro, Tema 3: Icke-linjära ekvationer, 18 november, 2015
(1 : 37)
Innehåll
I
Icke-linjära ekvationer och system
I
Anonyma funktioner Matlab
I
Konvergenshastighet
I
Newtons metod
I
Fixpunktsiteration
I
Minstakvadratproblem
Eddie Wadbro, Tema 3: Icke-linjära ekvationer, 18 november, 2015
(2 : 37)
Icke-linjära ekvationer
I
Många beräkningsproblem är till sin natur icke-linjära
I
Ex.: Eulers ekvationer i gasdynamiken (det aerodynamiska problem
som diskuterades i inledningen av första föreläsningen) blir, efter
diskretisering, ett mycket stort system av icke-linjära ekvationer som
beskriver tryck, densitet och hastigheter i noderna
I
Linjära ekvationssystem kan i princip lösas för hand. För stora
problem behövs datorer
I
Icke-linjära ekvationer kan sällan lösas exakt för hand. Det finns
ingen “formel” för lösningen i allmänhet. Det behövs numeriska
metoder även i det skalära fallet.
Eddie Wadbro, Tema 3: Icke-linjära ekvationer, 18 november, 2015
(3 : 37)
Icke-linjära system, exempel
x2
Exempel, sid. 94 i kursboken:
4x12 + 9x22 − 36 = 0
16x12 − 9x22 − 36 = 0
ska lösas för x1 ≥ 0, x2 ≥ 0
x1
Ovanligt exempel: kan lösas för hand!
Observera att problemet är linjärt i y1 = x12 , y2 = x22
Eddie Wadbro, Tema 3: Icke-linjära ekvationer, 18 november, 2015
(4 : 37)
Icke-linjära system
Ett system med n ekvationer och n okända skrivs generellt
f1 (x1 , . . . , xn ) = 0
f2 (x1 , . . . , xn ) = 0
..
.
fn (x1 , . . . , xn ) = 0
System ovan kan skrivas i vektorform som
f(x) = 0,
där
 
x1
 .. 
x=.
 
f1
 .. 
f =.
fn
och
xn
Eddie Wadbro, Tema 3: Icke-linjära ekvationer, 18 november, 2015
(5 : 37)
Icke-linjära system
För icke-linjära system så är Jacobimatrisen eller Jacobianen J av
central betydelse.
Jacobianens komponenter är
Jij =
∂fi
∂xj
vilket ger
 ∂f
1
∂x
 ∂f 1
 2
 ∂x1
J= .
 ..

∂fn
∂x1
∂f1
∂x2
∂f2
∂x2
∂fn
∂x2
...
...
...
∂f1
∂xn

∂f2 
∂xn 

∂fn
∂xn



För det inledande exemplet är Jacobianen
8x1
18x2
J=
32x1 −18x2
Eddie Wadbro, Tema 3: Icke-linjära ekvationer, 18 november, 2015
(6 : 37)
Linjär kontra icke-linjär
I
Linjära system har entydig lösning om systemmatrisen är inverterbar
I
För icke-linjära ekvationer är det svårt att verifiera om de har en
lösning innan någon beräkning har utförts
I
Icke-linjära system har ofta mer än en lösning
f (x)
2
1
Figuren till vänster illustrerar att
funktionen f (x) = cos(3x)e −x − x
x har tre rötter (lösningar till f (x) = 0)
i intervallet [−2, 2]
0
−1
−2
−2
−1
0
1
2
Eddie Wadbro, Tema 3: Icke-linjära ekvationer, 18 november, 2015
(7 : 37)
Icke-linjära ekvationer
Ex.: Lös den icke-linjära ekvationen cos(3x)e −x − x = 0. Denna ekvation
kan inte lösas “analytiskt,” så vi försöker beräkna en numerisk lösning
För att hitta en lösning x med utgånspunkt 0.5 i Matlab kan man skriva
>> format long
>> f = @(x) cos(3*x).*exp(-x)-x
f =
@(x)cos(3*x).*exp(-x)-x
>> x0 = 0.5;
>> fzero(f,x0)
ans =
0.350198594439993
I
I Matlab kan man definiera funktioner antingen i m-filer eller genom
att använda anonyma funktioner som ovan
I
@(x) i ovanstående kod instruerar Matlab att f är en funktion (och
inte ett tal eller en matris) och att x är den oberoende variabeln i
denna funktion
Eddie Wadbro, Tema 3: Icke-linjära ekvationer, 18 november, 2015
(8 : 37)
Anonyma funktioner i Matlab
Att skapa en anonym funktion är ett alternativ till att skriva och spara en
m-fil
Syntaxen för att skapa en anonym funktion är
fhandle = @(arglist) expr
Elementen i högerledet är
I
I
I
@ är Matlab operatorn som skapar funktionshandtaget
arglist är en kommaseparerad lista om innehåller inargumenten till
funktionen
expr representerar funktionskroppen, d.v.s., den kod som kommer
att köras när funktionen anropas
Observera. Funktionshandtag ger inte bara tillgång till anonyma funktioner. Genom att
använda en något modifierad syntax kan man skapa ett funktionshandtag till godtyckliga
Matlab funktioner. Ex.:
fhandle = @exp
Eddie Wadbro, Tema 3: Icke-linjära ekvationer, 18 november, 2015
(9 : 37)
Anonyma funktioner i Matlab
Anonyma funktioner innehåller vanligen variabler av två olika “sorter”
I
Variabler som specificeras i inargumentlistan
I
Variabler som anges i själva uttrycket. Matlab fångar dessa variables
när funktionen skapas och håller dem konstanta under hela
livslängden för funktionshandtaget
>> a=2;
>> f = @(x) cos(a*x).*exp(-x)-x;
>> f(1)
ans =
-1.153091865674226
>> a=3;
>> f(1)
ans =
-1.153091865674226
>> f = @(x) cos(a*x).*exp(-x)-x;
>> f(1)
ans =
-1.364197886413293
Sätt a till 2
Låt f vara ett funktionshandtag till
cos(ax)e −x − x, med nuvarande a, (a = 2)
evaluera f(1)
ändra a
evaluera f(1)
Låt f vara ett funktionshandtag till
cos(ax)e −x − x, med nuvarande a, (a = 3)
evaluera f(1)
Eddie Wadbro, Tema 3: Icke-linjära ekvationer, 18 november, 2015
(10 : 37)
Direkta kontra iterativa metoder
Numeriska metoder för ekvationslösning är antingen direkta eller
iterativa
I
Direkta metoder: den exakta lösningen (upp till avrudningfel) hittas
efter ett fixt (känt i förväg) antal aritmetiska operationer. Ex.:
Gausselimination för att lösa ett system av linjära ekvationer
I
Iterativa metoder: algoritmen producerar en sekvens av
approximationer x1 , x2 , . . . som (förhoppningsvis) närmar sig den
(okända) exakta lösningen x. Behöver avsluta algoritmen när
kxk − xk är tillräckligt liten; vilket introducerar ytterligare fel
I
För att lösa Icke-linjära ekvationer behövs generellt iterativa metoder
I
Iterativa metoder är nödvändiga även för att lösa mycket stora
linjära system, när direkta metoder är alltför tids och minneskrävande
Eddie Wadbro, Tema 3: Icke-linjära ekvationer, 18 november, 2015
(11 : 37)
Icke-linjära lösare i Matlab
I
fzero hittar ett nollställe till en kontinuerlig funktion av en variabel
I
Det finns ingen lösare för icke-linjära system i “grund” Matlab
I
fsolve löser icke-linjära ekvationssystem med hjälp av Newtons
metod, den mest använda metoden för att lösa små till medelstora
icke-linjära ekvationssystem. Denna funtion finns i Matlabs
optimization toolbox vilken säljs separat (ingår i universitetets licens
och finns tillgänglig i datasalarna)
Eddie Wadbro, Tema 3: Icke-linjära ekvationer, 18 november, 2015
(12 : 37)
Iterativa metoder
Tumregel:
icke-linjära ekvationer eller mycket stora linjära
ekvationssystem ⇒ iterativa metoder
En generell iterativ metod:
Vanligtvis:
x = startgissning
x 0 = startgissning
while stoppkriterium inte mött
while stoppkriterium inte mött
Hitta en riktning p k och
en steglängd αk
x = ny gissning
Sätt x k+1 = x k + αk p k
end
end
Eddie Wadbro, Tema 3: Icke-linjära ekvationer, 18 november, 2015
(13 : 37)
Konvergenshastighet
För linjära system så kan vi hitta den entydiga lösningen i ett ändligt
antal steg
Vanligtvis så producerar iterativa metoder enbart approximativa lösningar
⇒ Avsluta när lösningen är tillräckligt bra
Viktigt: Konvergenshastighet
Hur snabbt x k (approximationen vid steg k) närmar sig den exakta
lösningen x ∗
Eddie Wadbro, Tema 3: Icke-linjära ekvationer, 18 november, 2015
(14 : 37)
Konvergenshastighet
Definition
∗
Följden {x k }∞
k=1 konvergerar mot x med konvergensordning (Eng.
convergence rate) p och asymptotisk felkonstant (Eng. rate constant)
C < ∞ om x k → x ∗ och
kx k+1 − x ∗ k
lim
=C
k→∞ kx k − x ∗ kp
I
Linjär: p = 1 och 0 < C < 1. Felet multipliceras essentiellt med C
varje iteration
I
Kvadratisk: p = 2. Ungefär en fördubbling av antalet korrekta siffror
varje iteration
I
Superlinjär: p = 1 och C = 0. “Snabbare än linjär”. Inkluderar
kvadratisk konvergens men även“mellanliggande”
konvergensordningar
Observera: definitionen gäller asymptotisk konvergenshastighet
Eddie Wadbro, Tema 3: Icke-linjära ekvationer, 18 november, 2015
(15 : 37)
Iterativa metoder: konvergenshastighet
Ex.: Linjär konvergens
Antag att kx 0 − x ∗ k = 0.5 och
kx k+1 − x ∗ k = 0.3kx k − x ∗ k för varje k > 0
100 •
I
x ∗k
0.3n ,
kx n −
= 0.5 ·
alltså,
felet avtar exponentiellt
Fel
I
•
10−1
Då gäller,
Felet blir aldrig noll!
10
−2
10
−3
•
•
•
10−4
•
0
1
2
3
Iteration
Eddie Wadbro, Tema 3: Icke-linjära ekvationer, 18 november, 2015
4
5
(16 : 37)
Iterativa metoder: konvergenshastighet
Ex.: Kvadratisk konvergens
Antag att kx 0 − x ∗ k = 0.5 och
kx k+1 − x ∗ k = kx k − x ∗ k2 för varje k > 0
100 •
10
I
I
n
kx n − x ∗ k = 0.52 , alltså, antalet
korrekta siffror dubbleras i stort
sett varje iteration
Fel
Då gäller,
Felet blir aldrig noll!
•
•
−2
•
10−4
•
10−6
10−8
10−10
•
0
1
2
3
Iteration
4
5
Observera: Värdet på konstanten C spelar inte någon roll i definitionen av
kvadratisk konvergens, medan C < 1 är ett krav i definitionen av linjär
konvergens
Eddie Wadbro, Tema 3: Icke-linjära ekvationer, 18 november, 2015
(17 : 37)
Newtons metod för icke-linjära ekvationer
Vi vill lösa ekvationen f (x) = 0 där f är kontinuerligt deriverbar
Vidare, anta att vi är vid iteration k och att vår nuvarande gissning xk
inte löser vårt problem (f (xk ) 6= 0)
Vi vill hitta ett steg sk så att f (xk + sk ) ≈ 0
Genom att Taylor utveckla f kring xk får vi
f (xk + sk ) = f (xk ) + f 0 (xk )sk + O(|sk |2 )
Ide: välj sk = −f (xk )/f 0 (xk ) så att de två första termerna ovan tar ur
varandra. Sätt xk+1 = xk + sk
Eddie Wadbro, Tema 3: Icke-linjära ekvationer, 18 november, 2015
(18 : 37)
Newtons metod för en icke-linjär ekvation:
Alternativ tolkning
Taylor utveckling ger
f (xk + sk ) = f (xk ) + f 0 (xk )sk + O(|sk |2 )
Vilket ger att funktionen
fˆ(x) = f (xk ) + f 0 (xk )(x − xk )
approximerar f kring xk
Vi kan lösa den linjära ekvationen fˆ(x) = 0 för att hitta nästa gissning
xk+1
Eddie Wadbro, Tema 3: Icke-linjära ekvationer, 18 november, 2015
(19 : 37)
Newtons metod för icke-linjära ekvationer
n ekvationer fi (x) = 0, i = 1, . . . , n, i n okända x = (x1 , . . . , xn )
Vi antar att varje fi är kontinuerligt differentierbar
Antag att vid iteration k så har vi fi (xk ) 6= 0. Vi vill hitta ett steg
(k)
(k)
sk = (s1 , . . . , sn ) så att fi (xk + sk ) ≈ 0 för alla i
Taylor utveckling ger
fi (xk + s) = fi (xk ) +
n
X
∂fi
(xk )sj + . . .
∂xj
j=1
Ide: Välj s så att de två första termerna ovan tar ut varandra för varje i
fi (xk ) +
n
X
∂fi
(k)
(xk )sj = 0,
∂xj
i = 1, . . . , n
j=1
Eddie Wadbro, Tema 3: Icke-linjära ekvationer, 18 november, 2015
(20 : 37)
Newtons metod för icke-linjära ekvationer
Grundalgoritmen:
1. Välj en startgissning x0
2. For k = 0, 1, . . .
2.1 Lös det linjära ekvationssystemet
J(x k )s k = −f (x k ),
där J(x k ) är Jacobimatrisen för f evaluerad i punkten x k
2.2 Sätt
x k+1 = x k + s k
Newtons metod omvandlar problemet att lösa ett icke-linjärt system till
problemet att lösa en sekvens av linjära system
Eddie Wadbro, Tema 3: Icke-linjära ekvationer, 18 november, 2015
(21 : 37)
Newtons metod för icke-linjära ekvationer
I
Som ovan, låt f : Rn → Rn vara en kontinuerligt differentierbar
funktion
I
Antag att det finns en lösning till det icke-linjära systemet fi (x) = 0.
(f (x ∗ ) = 0 för något x ∗ ∈ Rn )
Sats
Följden x 0 , x 1 , . . . som genereras av Newtons metod konvergerar
kvadratiskt mot x ∗ om
1. J(x ∗ ) är inverterbar och
2. kx 0 − x ∗ k är tillräckligt liten
Eddie Wadbro, Tema 3: Icke-linjära ekvationer, 18 november, 2015
(22 : 37)
Egenskaper hos Newtons metod
I
Fördel: Kan konvergera extremt snabbt
I
Begänsningar:
1.
2.
3.
4.
Vi
Vi
Vi
Vi
behöver beräkna Jacobianen varje iteration
behöver lösa ett linjärt ekvationssystem varje iteration
behöver starta tillräckligt nära lösningen
kan få problem om Jacobianen är singular vid någon iteration
Eddie Wadbro, Tema 3: Icke-linjära ekvationer, 18 november, 2015
(23 : 37)
Jacobianberäkningar
I
En vanlig situation: funktionsvärdena beräknas med hjälp av någon
kommersiell mjukvara som inte ger tillgång till källkoden. Speciellt är
Jacobianen inte tillgänglig!
I
Jacobianen kan approximeras med finita differenser:
∂fi
(x1 , . . . , xn )
∂xj
fi (x1 , . . . , xj + h, . . . , xn ) − fi (x1 , . . . , xj , . . . , xn )
≈
h
Jij (x1 , . . . , xn ) =
I
h ska inte vara alltför stor (dålig approximation av derivaten) eller för
liten (kancellering av signifikanta siffror). Bästa värdet på h är
√
problemberoende. Tumregel: h ∼ M
Eddie Wadbro, Tema 3: Icke-linjära ekvationer, 18 november, 2015
(24 : 37)
Jacobianberäkningar
Beräkning av Jacobimatrisen med finita differenser kräver n2 + 1
evalueringar av f .
I
Kan vara mycket kostsamt för stora problem
I
(speciellt om en evaluering av f innefattar en numerisk lösning av
exempelvis en partiell differentialekvation)
Ett vanligt tillvägagånssätt
I
Börja med finita differensapproximationer av Jacobianen
I
Använd sedan sekantapproximationer, baserade på
funktionsvärdena f (x k+1 ) och f (x k ) samt gissningarna x k+1 och
x k , för att uppdatera Jacobianapproximationen (Begränsning 1).
I
Användning av sekantapproximationer istället för den exakta
Jacobianen reducerar konvergensordningen något (metoden
konvergerar superlinjärt istället för kvadratiskt)
Eddie Wadbro, Tema 3: Icke-linjära ekvationer, 18 november, 2015
(25 : 37)
Andra problem med Newtons metod
I
Behovet av en startgissning nära lösningen kan tas bort genom
användning av globaliseringstekniker, vilket medför att metoden
konvergerar oavsett startgissning.
I
I
Globaliseringen modifierar steglängden och/eller stegriktningen när
man är långt från lösningen, vilket är typiskt under de första stegen.
(Begränsning 3)
Det finns också tekniker för att ta hand om singulära Jacobianer
(Begränsning 4)
Spjutspetsimplementationer av Newtons metod är mycket mer
komplicerade än grundalgoritmen!
Eddie Wadbro, Tema 3: Icke-linjära ekvationer, 18 november, 2015
(26 : 37)
Fixpunktsiterationer
I
En klass av iterativa metoder som undviker lösning av linjära system
(Begränsning 2)
I
De är därför extra intressanta för lösning av mycket stora problem
I
Vanligtvis mycket långsammare än Newtons metod
Obervera: konstruktionen av en bra fixpunktsmetod är problemberoende
Eddie Wadbro, Tema 3: Icke-linjära ekvationer, 18 november, 2015
(27 : 37)
Fixpunktsiterationer
Ide:
I
Skriv f (x) = 0 som ett fixpunktsproblem x = g (x)
I
I
Det går alltid att göra en sådan omskrivning, exempelvis x = x − f (x)
Definiera det iterativa schemat
x n+1 = g (x n )
och “hoppas” att det konvergerar mot en fixpunkt till g , vilket per
konstruktion är en lösning till f (x) = 0
Eddie Wadbro, Tema 3: Icke-linjära ekvationer, 18 november, 2015
(28 : 37)
Fixpunktsiterationer
I
Iterationerna konvergerar om g är en kontraktion; alltså om det
finns ett C < 1 så att
kg (x) − g (y )k ≤ C kx − y k
för varje x, y i en (konvex) delmängd av Rn
I
Eller i ord, om g (x) och g (y ) är strikt närmare varandra än x och y
så är g en kontraktion
Eddie Wadbro, Tema 3: Icke-linjära ekvationer, 18 november, 2015
(29 : 37)
Fixpunktsiterationer
I
Inte alltid enkelt att se om en avbildning är en kontraktion!
I
Ett tillräckligt villkor för att en kontinuerligt differentierbar funktion
g ska vara en kontraktion är att
kJ(x)k < 1
för varje x i ett område, där J är Jacobianen för g (Observera:
matrisnorm)
Eddie Wadbro, Tema 3: Icke-linjära ekvationer, 18 november, 2015
(30 : 37)
Fixpunktsiterationer
I
Enkelt, vi behöver inte lösa några linjära ekvationssystem
I
För extremt stora problem kan de vara den enda möjliga
lösningsvägen
I
Att konstruera en lämplig kontraktiv funktion g är i hög grad
problemberoende
I
Robusta versioner av Newtons metod är (om tillämpliga) vanligtvis
mycket snabbare
Eddie Wadbro, Tema 3: Icke-linjära ekvationer, 18 november, 2015
∗
(31 : 37)
Minstakvadratproblem
Antag att vi har observationer
yi
i = 1, 2, . . . . , m vid tiderna
ti
i = 1, 2, . . . . , m
◦
och en model y (t, x)
◦ ◦
◦ ◦
◦
Ex1: y (t, x) = x1 + x2 t + x3 t
2
Ex2: y (t, x) = x1 e x2 t
Vi vill bestämma de “bästa” x1 , x2 , . . . , xn så att y (ti , x) ≈ yi
(Vanligtvis m >> n)
Eddie Wadbro, Tema 3: Icke-linjära ekvationer, 18 november, 2015
∗
(32 : 37)
Minstakvadratproblem
Vanligtvis (på grund av mätfel, modelleringsfel, brus, etcetera),
yi − y (ti , x) 6= 0 även för det “bästa” parametervalet x
Här betyder “bästa” parametervalet x det x som är bäst i
minstakvadratmening
Vi vill alltså:
Hitta det x ∈ Rn som minimerar
M
1X
[yi − y (ti , x)]2
2
i=1
Definition
Linjärt minstakvadratproblem: y är linjär i x (Ex1)
Icke-linjärt minstakvadratproblem: y är icke-linjär i x (Ex2)
Eddie Wadbro, Tema 3: Icke-linjära ekvationer, 18 november, 2015
(33 : 37)
∗
Linjära minstakvadratproblem
Låt oss först titta på fallet
m
1X
minn f (x) =
fi (x)2
x∈R
2
i=1
där fi är linjär i x, vilket betyder att
2 1
1
1
1
f (x) = Ax−b = (Ax−b)(Ax−b)T = x T AT Ax−x T AT b+ b T b,
2
2
2
2
där A är en m × n matris och b är en m kolonnvektor
Om AT A är positivt definit så har minimeringsproblemet ovan lösningen
x ∗ , som uppfyller normalekvationerna
AT Ax ∗ = AT b
Eddie Wadbro, Tema 3: Icke-linjära ekvationer, 18 november, 2015
∗
(34 : 37)
Icke-linjära minstakvadratproblem
m
minn f (x) =
x∈R
1X
fi (x)2
2
i=1
I detta fall kan vi skriva gradienten och Hessianen av f som
∇f (x) =
∇2 f (x) =
m
X
∇fi (x)fi (x)
i=1
m h
X
∇fi (x) ∇fi (x)
T
i
+ ∇2 fi (x)fi (x)
i=1
Vi kan exempelvis använda oss av Newtons metod för att lösa ∇f = 0
Eddie Wadbro, Tema 3: Icke-linjära ekvationer, 18 november, 2015
∗
(35 : 37)
Icke-linjära minstakvadratproblem
m
minn f (x) =
x∈R
1X
fi (x)2
2
i=1
Det är dock vanligtvis bättre att använda den struktur som problemet har
Taylor utveckling ger
1
f (x + s) = f (x) + ∇f (x)T s + s T ∇2 f (ξ)s
2
där ξ = x + αs för något α ∈ [0, 1]
Ide: vid iteration k, approximera Hessianen ∇2 f i den okända punkten ξ
m h
med
i
X
T
Bk =
∇fi (x k ) ∇fi (x k ) + någonting enkelt
i=1
Eddie Wadbro, Tema 3: Icke-linjära ekvationer, 18 november, 2015
(36 : 37)
∗
Icke-linjära minstakvadratproblem
Vad är “någonting enkelt”?
I
0 (noll) rimligt om redidualen är liten samt för svagt olinjära problem
I
(annars använd λI , där I är identitetsmatrisen och λ > 0)
Approximera f med
1
fˆ(x) = f (x k ) + ∇f (x k )T (x − x k ) + (x − x k )T Bk (x − x k )
2
ta fram steget s k genom att lösa Bk s k = −∇f
Observera:
Ekvivalent till att lösa det linjära minstakvadratproblemet
∇F (∇F )T s = −(∇F )F
där F = [f1 , f2 , . . . , fn ] och ∇F = [∇f1 ∇f2 . . . ∇fn ]T
Eddie Wadbro, Tema 3: Icke-linjära ekvationer, 18 november, 2015
(37 : 37)