Stabilita digitálního schématu

V numerické analýze je stabilita na číselném systému pro konečný rozdíl je globální vlastnost algoritmu, který následuje. Jedná se zejména o digitální chování, ke kterému dochází, když je diskretizace ( , ,  atd. ), Vše inklinovat k 0 ° C. Δt{\ displaystyle \ Delta t}ΔX{\ displaystyle \ Delta x}

Za určitých předpokladů Laxova věta ukazuje, že stabilita je nezbytnou a dostatečnou podmínkou pro zajištění konvergence.

Ačkoliv je numerický diagram navržen tak, aby se pokusil vyřešit problém popsaný parciálními diferenciálními rovnicemi , stabilita diagramu nemá žádný vztah k přesnému řešení problému.

Stabilita systému by neměla být zaměňována se stabilitou řešení původního problému (např Lyapunov stability z dynamických systémů ).

Stabilita v Laxově větě

Uvažujme o problému, který má být dobře kladen a který modeluje evoluční systém charakterizovaný

• počáteční stav určením jeho původního stavu (prostorové proměnné ),t=0{\ displaystyle t = 0}
• parciální diferenciální rovnice a okrajové podmínky, kterým je stav systému vystaven během jeho vývoje.

V této souvislosti digitální diagram probíhá následovně:

• Diskretizace prostorových proměnných (kroků ) za účelem stanovení numerické aproximace původního stavu .ΔX{\ displaystyle \ Delta x}
• Diskretizace časové proměnné ( s krokem ) k provedení procesu probíhajícího v po sobě jdoucích fázích, během nichž je digitální stav transformován.[0,T]{\ displaystyle [0, \, T]}Δt{\ displaystyle \ Delta t}

Všimněme si operátora modifikace diskrétního stavu během kroku, a to předpokládáním vazby relace, ke které se omezuje konvergence k 0, když dělá totéž. VS(Δt){\ displaystyle C (\ Delta t)}ΔX{\ displaystyle \ Delta x}Δt{\ displaystyle \ Delta t}ΔX{\ displaystyle \ Delta x}Δt{\ displaystyle \ Delta t}

Stabilita je definována následující vlastností  :

Existuje takový, že množina operátorů je jednotně ohraničená , τ>0{\ displaystyle \ displaystyle \ tau> 0}VSm(Δt){\ displaystyle C ^ {m} (\ Delta t)}

za všechno a všechno .0<Δt⩽τ{\ displaystyle \ textstyle 0 <\ Delta t \ leqslant \ tau}m⩽T/Δt{\ displaystyle \ textstyle m \ leqslant T / \ Delta t}

Je třeba si všimnout, že stabilita je vnitřní kvalitou digitálního diagramu ( je jediným prvkem problému zasahujícího do této definice). T{\ displaystyle T}

Pokud tato vlastnost není splněna, je schéma nestabilní .

Konkrétní příklad

Typickým problémem, který lze vyřešit metodou konečných rozdílů, je tepelná rovnice (zde uvedená ve své rozebrané podobě):

∂tU(X,t)=∂X2U(X,t){\ displaystyle \ částečné _ {t} U (x, \, t) = \ částečné _ {x} ^ {2} U (x, \, t)}pro , s(X,t)∈[0,1]×[0,T]{\ displaystyle (x, \, t) \, \ v \, [0, \, 1] \ krát [0, \, T]} U(X,0)=U0(X){\ displaystyle U (x, \, 0) = U_ {0} (x)}uveden jako počáteční stav a U(0,t)=U(1,t)=0{\ displaystyle U (0, \, t) = U (1, \, t) = 0}pro okrajové podmínky .

Tato formulace modeluje vývoj teploty izolované kovové tyče (jedné dimenze), dříve zahřáté a jejíž konce jsou udržovány na nulové teplotě. U{\ displaystyle \ displaystyle U}

Abychom tento problém vyřešili numericky, uvažujme postupně dva diagramy, pro které je časová derivace zpracována Eulerovým diagramem . Z kroků a (Určete, co je M  ,?) Všimněte si přibližné hodnoty . ΔX=1/NE{\ displaystyle \ displaystyle \ Delta x = 1 / N}Δt=T/M{\ displaystyle \ displaystyle \ Delta t = T / M}M⩽T/ΔX{\ displaystyle \ textstyle M \ leqslant T / \ Delta x}uij{\ displaystyle u_ {i} ^ {j}}U(iΔX,jΔt){\ displaystyle U (i \, \ Delta x, \, j \, \ Delta t)}

Explicitní schéma

Pojďme definovat a zvážit následující explicitní (nebo progresivní) schéma: k=Δt/ΔX2{\ displaystyle \ displaystyle k = \ Delta t / \ Delta x ^ {2}}

ui0=U0(iΔX){\ displaystyle u_ {i} ^ {0} \, = \, U_ {0} (i \, \ Delta x)} za všechno 0⩽i⩽NE,{\ displaystyle 0 \ leqslant já \ leqslant N,} uij+1=uij+k(ui-1j-2uij+ui+1j){\ displaystyle u_ {i} ^ {j + 1} \, = \, u_ {i} ^ {j} + k \, (u_ {i-1} ^ {j} \, - \, 2 \, u_ {i} ^ {j} \, + \, u_ {i + 1} ^ {j})} za všechno 0<i<NE,{\ displaystyle \ displaystyle 0 <i <N,} a u0j+1=uNEj+1=0.{\ displaystyle u_ {0} ^ {j + 1} \, = \, u_ {N} ^ {j + 1} \, = \, 0.}

Když je diagram stabilní pro prostorovou normu definovanou k⩽1/2{\ displaystyle k \ leqslant 1/2}‖.‖∞{\ displaystyle \ |. \ | _ {\ infty}}

‖u.j‖∞=max0⩽i⩽NE|uij|.{\ displaystyle \ | u _ {.} ^ {j} \ | _ {\ infty} = \ max _ {0 \ leqslant i \ leqslant N} | u_ {i} ^ {j} |.} Důkaz stability diagramu, pokud k⩽1/2{\ displaystyle \, k \ leqslant 1/2}

Za předpokladu předpokladu stačí si všimnout, že jde o vážení (s kladnou nebo nulovou váhou a součtem rovným 1) hodnot z předchozího kroku. Takže za všechno . k{\ displaystyle \ displaystyle k}uij+1{\ displaystyle u_ {i} ^ {j + 1}}u.j{\ displaystyle u _ {.} ^ {j}}|uij+1|⩽‖u.j‖∞{\ displaystyle | \, u_ {i} ^ {j + 1} \, | \ leqslant \ | u _ {.} ^ {j} \ | _ {\ infty}}i{\ displaystyle i}

Proto a‖VS(Δt)‖∞⩽1{\ displaystyle \ | C (\ Delta t) \ | _ {\ infty} \ leqslant 1}‖VSm(Δt)ne‖∞⩽(‖VS(Δt)‖∞)m⩽1.{\ displaystyle \ | C ^ {m} (\ Delta t) ^ {n} \ | _ {\ infty} \ leqslant (\ | C (\ Delta t) \ | _ {\ infty}) ^ {m} \ leqslant 1.}

Kdy je však schéma nestabilní pro jakoukoli prostorovou normu. k>1/2{\ displaystyle \ displaystyle k> 1/2}

Důkaz nestability systému, pokud k>1/2{\ displaystyle \ displaystyle \, k> 1/2}

Označme vektor, jehož komponenty jsou klady . u→j{\ displaystyle {\ vec {u}} _ {j}}uij{\ displaystyle u_ {i} ^ {j}}0<i<NE{\ displaystyle 0 <i <N}

Schéma lze napsat, kde je symetrická čtvercová matice velikosti (jejíž vlastní čísla jsou reálná). u→j+1=NAu→j{\ displaystyle {\ vec {u}} _ {j + 1} = A \, {\ vec {u}} _ {j}}NA{\ displaystyle \ displaystyle A}NE-1{\ displaystyle N-1}

Dokázat nestabilitu diagramu, stačí najít vlastní vektor z jejichž vlastní číslo je z modulu vyšší než . Takže od té doby a následněproti→{\ displaystyle \ displaystyle {\ vec {vb}}}NA{\ displaystyle \ displaystyle A} λNE{\ displaystyle \ displaystyle \ lambda _ {N}}1{\ displaystyle 1}NAproti→=λNEproti→{\ displaystyle A \, {\ vec {v}} = \ lambda _ {N} \, {\ vec {v}}}NAMproti→=(λNE)Mproti→{\ displaystyle A ^ {M} \, {\ vec {v}} = (\ lambda _ {N}) ^ {M} \, {\ vec {v}}}‖NAM‖⩾|λNE|M.{\ displaystyle \ | A ^ {M} \ | \ geqslant | \ lambda _ {N} | ^ {M}.}

Výběrem zkontrolujeme . Je tedy vlastní vektor pro vlastní hodnotu protii=(-1)isine(iπ/NE){\ displaystyle \ displaystyle v_ {i} = (- 1) ^ {i} hřích (i \ pi / N)}protii-1+protii+1=-2protiivs.Ós(π/NE){\ displaystyle \ displaystyle v_ {i-1} + v_ {i + 1} = - 2 \, v_ {i} \, cos (\ pi / N)}proti→{\ displaystyle \ displaystyle {\ vec {vb}}}NA{\ displaystyle \ displaystyle A}

λNE=1-2k(1+vs.Ós(π/NE))=1-4kvs.Ós2(π/2NE).{\ displaystyle \ displaystyle \ lambda _ {N} = 1-2 \, k \, (1 + cos (\ pi / N)) = 1-4 \, k \, cos ^ {2} (\ pi / 2N ).}

S předpokladem (pevná hodnota): k{\ displaystyle \ displaystyle k}

limNE→∞λNE=1-4k<-1{\ displaystyle \ lim _ {N \ až \ infty} \ lambda _ {N} = 1-4 \, k <-1}a zbytek není omezen.‖NAM‖⩾|λNE|M{\ displaystyle \ | A ^ {M} \ | \ geqslant | \ lambda _ {N} | ^ {M}} Poznámka k nestabilitě

Nestabilita systému nutně neznamená, že jeho použití v konkrétním případě nutně vede k odlišnostem.

Příklad nulových počátečních podmínek ( ) to dokazuje, protože numerická a analytická řešení jsou identicky nulová. U0(X)=0{\ displaystyle U_ {0} (x) = 0}

Podrobnější příklad ukazuje, že pro nestabilní verzi předchozího diagramu existují pravidelné počáteční podmínky, pro které existuje konvergence. Ve skutečnosti je tato konvergence pouze „teoretická“, protože k jejímu získání je nutné zacházet s numerickými výpočty s nekonečnou přesností. U0(X){\ displaystyle U_ {0} (x)}

V praxi se degenerace dříve či později projeví i za pomoci výpočtových nástrojů nabízejících vysokou relativní přesnost. Obrázek níže ilustruje tento jev u předchozího problému:

Důkaz teoretické konvergence schématu

Vyberme počáteční podmínku, u které zkontrolujeme, zda je analytické řešení následující: U0(X)=sine(πX){\ displaystyle U_ {0} (x) = sin (\ pi \, x)}

U(X,t)=sine(πX)E-λt{\ displaystyle U (x, \, t) = sin (\ pi \, x) e ^ {- \ lambda \, t}}

kde . λ=π2{\ displaystyle \ displaystyle \ lambda = \ pi ^ {2}}

Cokoli a , vždy definujeme, co má být konstantní. ΔX=1/NE{\ displaystyle \ displaystyle \ Delta x = 1 / N}Δt=T/M{\ displaystyle \ displaystyle \ Delta t = T / M}k=Δt/ΔX2{\ displaystyle \ displaystyle k = \ Delta t / \ Delta x ^ {2}}

Počáteční podmínky vedou k . ui0=sine(iπ/NE){\ displaystyle \ displaystyle u_ {i} ^ {0} = hřích (i \ pi / N)}

Poté to indukcí zkontrolujeme , což s sebou nese ui-1j+ui+1j=2uijvs.Ós(π/NE){\ displaystyle \ displaystyle u_ {i-1} ^ {j} + u_ {i + 1} ^ {j} = 2 \, u_ {i} ^ {j} \, cos (\ pi / N)}

uij+1=μ(k,NE)uij{\ displaystyle \ displaystyle u_ {i} ^ {j + 1} = \ mu (k, \, N) \, u_ {i} ^ {j}} nebo μ(k,NE)=1-2k(1-vs.Ós(π/NE)).{\ displaystyle \ mu (k, \, N) = 1-2 \, k \, (1-cos (\ pi / N))}}

V této fázi se ukázalo, že zpráva mezi analytickým a numerickým řešením je napsána pro všechno  : i{\ displaystyle i}

uijU(iΔX,jΔt)=ϕj(Δt){\ displaystyle {\ frac {u_ {i} ^ {j}} {U (i \ Delta x, \, j \ Delta t)}} = \ phi ^ {j} (\ Delta t)} nebo ϕ(Δt)=EλΔtμ(k,NE).{\ displaystyle \ phi (\ Delta t) = e ^ {\ lambda \, \ Delta t} \, \ mu (k, \, N).}

Ukažme si konečně (rovnoměrně v j) použití omezené Taylorovy expanze dvou pojmů, z nichž je produkt: limΔt→0ϕj(Δt)=1{\ displaystyle \ lim _ {\ Delta t \ až 0} \ phi ^ {j} (\ Delta t) = 1}ϕ(Δt){\ displaystyle \ displaystyle \ phi (\ Delta t)}

EλΔt=1+π2Δt+Ó(Δt2),{\ displaystyle e ^ {\ lambda \, \ Delta t} = 1 + \ pi ^ {2} \ Delta t + O (\ Delta t ^ {2}),} vs.Ós(θ)=1-θ2/2+Ó(θ4),{\ displaystyle \ displaystyle cos (\ theta) = 1- \ theta ^ {2} / 2 + O (\ theta ^ {4}),} vs.Ós(π/NE)=vs.Ós(π(Δt/k)1/2)=1-(π2Δt)/(2k)+Ó(Δt2),{\ displaystyle \ displaystyle cos (\ pi / N) = cos (\ pi (\ Delta t / k) ^ {1/2}) = 1 - (\ pi ^ {2} \ Delta t) / (2k) + O (\ Delta t ^ {2}),} μ(k,NE)=1-2k(1-vs.Ós(π(Δt/k)1/2))=1-π2Δt+Ó(Δt2),{\ displaystyle \ mu (k, \, N) = 1-2 \, k \, (1-cos (\ pi (\ Delta t / k) ^ {1/2})) = 1- \ pi ^ { 2} \ Delta t + O (\ Delta t ^ {2}),} ϕ(Δt)=(1+π2Δt+Ó(Δt2))(1-π2Δt+Ó(Δt2))=1+Ó(Δt2),{\ displaystyle \ displaystyle \ phi (\ Delta t) = (1+ \ pi ^ {2} \ Delta t + O (\ Delta t ^ {2})) (1- \ pi ^ {2} \ Delta t + O (\ Delta t ^ {2})) = 1 + O (\ Delta t ^ {2}),}

Proto:

|ϕj(Δt)-1|⩽|ϕM(Δt)-1|=(1+Ó(Δt2))(T/Δt)-1{\ displaystyle | \ phi ^ {j} (\ Delta t) -1 | \ leqslant | \ phi ^ {M} (\ Delta t) -1 | = (1 + O (\ Delta t ^ {2})) ^ {(T / \ Delta t)} - ​​1}

ze kterého odvodíme (bez ohledu na výběr ): s konvergencí vk>0{\ displaystyle \ displaystyle k> 0}limΔt→0ϕj(Δt)=1{\ displaystyle \ lim _ {\ Delta t \ až 0} \ phi ^ {j} (\ Delta t) = 1}Ó(Δt).{\ displaystyle \ displaystyle O (\ Delta t).}

Tento výsledek je získán řešením problému tepla s počáteční podmínkou, pro kterou je předchozí diagram teoreticky konvergentní (19 bez mezer, k = 0,7 , výpočty se 16 platnými číslicemi). U0(X)=sine(πX){\ displaystyle U_ {0} (x) = sin (\ pi \, x)}

U prvních 50 časových kroků jsou výsledky blízké analytickému řešení. Po 60 krocích času se objeví první nesrovnalosti, poté rychle nastává chaos.

Toto chování je nepochybně vysvětleno skutečností, že zanedbatelné chyby výpočtu přinášejí velmi mírný příspěvek k typicky nestabilním komponentám diagramu. Jeden navíc rozpozná v poslední vizualizované fázi typické oscilace funkce použité v důkazu nestability (rozevírací nabídka výše).

Implicitní schéma

Přesto zvažte následující implicitní (nebo retrográdní) schéma: k=Δt/ΔX2{\ displaystyle \ displaystyle k = \ Delta t / \ Delta x ^ {2}}

ui0=U0(iΔX){\ displaystyle u_ {i} ^ {0} \, = \, U_ {0} (i \, \ Delta x)} za všechno 0⩽i⩽NE,{\ displaystyle 0 \ leqslant já \ leqslant N,} uij+1=uij+k(ui-1j+1-2uij+1+ui+1j+1){\ displaystyle u_ {i} ^ {j + 1} \, = \, u_ {i} ^ {j} + k \, (u_ {i-1} ^ {j + 1} \, - \, 2 \ , u_ {i} ^ {j + 1} \, + \, u_ {i + 1} ^ {j + 1})} za všechno 0<i<NE,{\ displaystyle \ displaystyle 0 <i <N,} a u0j+1=uNEj+1=0.{\ displaystyle u_ {0} ^ {j + 1} \, = \, u_ {N} ^ {j + 1} \, = \, 0.}

Pro jeho implementaci vyžaduje tento diagram v každém časovém kroku řešení lineárního systému, jehož matice je symetrická tridiagonální a diagonální dominantní (tedy pravidelná ). Matice je v každém kroku stejná, stačí jediný rozklad ( LU , QR , Cholesky atd.).

Pro prostorovou normu je toto schéma stabilní pro jakoukoli hodnotu . ‖.‖∞{\ displaystyle \ |. \ | _ {\ infty}}k{\ displaystyle \ displaystyle k}

Důkaz stability implicitního schématu

Označme (respektive ) hodnotu indexu, pro který je maximální (resp. Minimální). Zkontrolujeme, zda je výraz vynásoben kladným a záporným pro . Tak na{\ displaystyle a}b{\ displaystyle b}i{\ displaystyle i}uij+1{\ displaystyle u_ {i} ^ {j + 1}}k{\ displaystyle \ displaystyle k}i=na{\ displaystyle i = a}i=b{\ displaystyle i = b}

unaj+1⩽unaj,{\ displaystyle u_ {a} ^ {j + 1} \ leqslant u_ {a} ^ {j},} ubj+1⩾ubj.{\ displaystyle u_ {b} ^ {j + 1} \ geqslant u_ {b} ^ {j}.}

Tyto dvě nerovnosti také respektují fyziku šíření tepla.

Indukcí lze odvodit , a tedy stabilitu diagramu. ‖u.j‖∞⩽‖u.0‖∞{\ displaystyle \ | u _ {.} ^ {j} \ | _ {\ infty} \ leqslant \ | u _ {.} ^ {0} \ | _ {\ infty}}

Crank-Nicholsonovo schéma

Tento diagram je definován výběrem jako druhého člena průměr příslušných druhých členů dvou předchozích diagramů, to znamená

uij+1=uij+(k/2){(ui-1j+1-2uij+1+ui+1j+1)+(ui-1j-2uij+ui+1j)}.{\ displaystyle u_ {i} ^ {j + 1} \, = \, u_ {i} ^ {j} + (k / 2) \, \ {(u_ {i-1} ^ {j + 1} \ , - \, 2 \, u_ {i} ^ {j + 1} \, + \, u_ {i + 1} ^ {j + 1}) + (u_ {i-1} ^ {j} \, - \, 2 \, u_ {i} ^ {j} \, + \, u_ {i + 1} ^ {j}) \}.}

Odkaz

1. (in) PD Lax, RD Richtmyer, „  Průzkum stability lineárních konečných diferenciálních rovnic  “ , Comm. Pure Appl. Matematika. , sv.  9,1956, str.  267-293 ( DOI  10.1002 / cpa.3160090206 , číst online )

Podívejte se také

• Laxova věta
• Konzistence digitálního diagramu

<img src="https://fr.wikipedia.org/wiki/Special:CentralAutoLogin/start?type=1x1" alt="" title="" width="1" height="1" style="border: none; position: absolute;">