Tuhá rovnice - Stiff equation

v matematika, a tuhá rovnice je diferenciální rovnice pro které jisté numerické metody pro řešení rovnice jsou numericky nestabilní, pokud není velikost kroku považována za extrémně malou. Ukázalo se, že je obtížné formulovat přesnou definici tuhosti, ale hlavní myšlenkou je, že rovnice obsahuje některé výrazy, které mohou vést k rychlým změnám v řešení.

Při numerické integraci diferenciální rovnice by se dalo očekávat, že požadovaná velikost kroku bude v oblasti, kde křivka řešení zobrazuje mnoho variací a je relativně velký tam, kde se křivka řešení narovnává a přibližuje se k přímce se sklonem téměř nulovým. U některých problémů tomu tak není. Aby numerická metoda poskytla spolehlivé řešení diferenciálního systému, je někdy nutné, aby velikost kroku byla na nepřijatelně malé úrovni v oblasti, kde je křivka řešení velmi hladká. Tento jev je znám jako ztuhlost. V některých případech mohou existovat dva různé problémy se stejným řešením, ale jeden není tuhý a druhý ano. Fenomén tedy nemůže být vlastností přesného řešení, protože je pro oba problémy stejný, a musí být vlastnictvím samotného diferenciálního systému. Takové systémy jsou tedy známé jako tuhé systémy.

Motivující příklad

Explicitní numerické metody vykazující nestabilitu při integraci tuhé obyčejné diferenciální rovnice

Zvažte problém počáteční hodnoty

${ Displaystyle , y '(t) = - 15 let (t), quad t geq 0, y (0) = 1.}$

(1)

Přesné řešení (zobrazené v azurové barvě) je

${ displaystyle y (t) = e ^ {- 15t} ,}$ s ${ displaystyle y (t) až 0}$ tak jako ${ displaystyle t to infty.}$

(2)

Hledáme a numerické řešení který vykazuje stejné chování.

Obrázek (vpravo) ilustruje numerické problémy pro různé numerické integrátory aplikované na rovnici.

Jeden z nejvýznamnějších příkladů tuhého Obyčejné diferenciální rovnice (ODR) je systém, který popisuje chemická reakce Robertsona^[1]:

${ displaystyle { dot {x}} = - 0,04x + 10 ^ {4} y cdot z}$ ${ displaystyle { dot {y}} = 0,04x-10 ^ {4} y cdot z-3 cdot 10 ^ {7} y ​​^ {2}}$ ${ displaystyle { dot {z}} = 3 cdot 10 ^ {7} y ​​^ {2}}$

(4)

Pokud někdo zachází s tímto systémem v krátkém intervalu, například ${ displaystyle t v [0,40]}$ v numerické integraci není problém. Pokud je však interval velmi velký (10¹¹ řekněme), pak mnoho standardních kódů nedokáže správně integrovat.

Dalšími příklady jsou sady ODR vyplývající z časové integrace velkých mechanismů chemických reakcí. Zde vzniká tuhost z koexistence velmi pomalých a velmi rychlých reakcí.^{[Citace je zapotřebí ]} K jejich vyřešení jsou softwarové balíčky KPP a Autochem může být použito.

Poměr tuhosti

Zvažte lineární konstantní koeficient nehomogenní systém

${ displaystyle mathbf {y} '= mathbf {A} mathbf {y} + mathbf {f} (x),}$

(5)

kde ${ displaystyle mathbf {y}, mathbf {f} in mathbb {R} ^ {n}}$ a ${ displaystyle mathbf {A}}$ je konstantní, diagonalizovatelný, ${ displaystyle n krát n}$ matice s vlastními hodnotami ${ displaystyle lambda _ {t} in mathbb {C}, t = 1,2, ldots, n}$ (předpokládá se odlišné) a odpovídající vlastní vektory ${ displaystyle mathbf {c} _ {t} v mathbb {C} ^ {n}, t = 1,2, ldots, n}$. Obecné řešení (5) má formu

${ displaystyle mathbf {y} (x) = součet _ {t = 1} ^ {n} kappa _ {t} exp ( lambda _ {t} x) mathbf {c} _ {t} + mathbf {g} (x),}$

(6)

kde κ_t jsou libovolné konstanty a ${ displaystyle mathbf {g} (x)}$ je zvláštní integrál. Nyní předpokládejme, že

${ displaystyle Re ( lambda _ {t}) <0, qquad t = 1,2, ldots, n,}$

(7)

což znamená, že každý z těchto výrazů${ displaystyle exp ( lambda _ {t} x) mathbf {c} _ {t} rightarrow 0}$tak jako ${ displaystyle x rightarrow infty}$, takže řešení ${ displaystyle mathbf {y} (x)}$přístupy ${ displaystyle mathbf {g} (x)}$ asymptoticky jako ${ displaystyle x rightarrow infty}$; termín ${ displaystyle exp ( lambda _ {t} x) mathbf {c} _ {t}}$ se bude monotónně rozpadat, pokud λ_t je reálné a sinusově, pokud λ_t je složité. tlumočení X být čas (jak to často bývá u fyzických problémů), ${ displaystyle Sigma _ {t = 1} ^ {n} kappa _ {t} exp ( lambda _ {t} x) mathbf {c} _ {t}}$ se nazývá přechodné řešení a ${ displaystyle mathbf {g} (x)}$ the ustálené řešení.Li ${ displaystyle | Re ( lambda _ {t}) |}$ je velký, pak odpovídající termín ${ displaystyle kappa _ {t} exp ( lambda _ {t} x) mathbf {c} _ {t}}$ rychle se rozpadne jakoX se zvyšuje a nazývá se tedy a rychlé přechodné; -li${ displaystyle | Re ( lambda _ {t}) |}$ je malý, odpovídající výraz${ displaystyle kappa _ {t} exp ( lambda _ {t} x) mathbf {c} _ {t}}$ rozpadá se pomalu a nazývá se a pomalý přechod. Nechat ${ displaystyle { overline { lambda}}, { podtržení { lambda}} in { lambda _ {t}, t = 1,2, ldots, n }}$ být definován

${ displaystyle | Re ({ overline { lambda}}) | geq | Re ( lambda _ {t}) | geq | Re ({ underline { lambda}}) |, qquad t = 1 , 2, ldots, n}$

(8)

aby ${ displaystyle kappa _ {t} exp ({ overline { lambda}} x) mathbf {c} _ {t}}$ je nejrychlejší přechodný a ${ displaystyle kappa _ {t} exp ({ podtržení { lambda}} x) mathbf {c} _ {t}}$ nejpomalejší. Nyní definujeme poměr tuhosti tak jako

${ displaystyle { frac {| Re ({ overline { lambda}}) |} {| Re ({ underline { lambda}}) |}}.}$[2]

(9)

Charakterizace tuhosti

V této části uvažujeme o různých aspektech fenoménu tuhosti. „Fenomén“ je pravděpodobně vhodnější slovo než „vlastnost“, protože to spíše naznačuje, že tuhost lze definovat přesnými matematickými termíny; ukázalo se, že to není možné uspokojivým způsobem, dokonce ani pro omezenou třídu systémů lineárních konstantních koeficientů. Uvidíme také několik kvalitativních výroků, které lze (a většinou byly) učiněny ve snaze zapouzdřit pojem tuhosti a co je pravděpodobně nejuspokojivější z nich, uvedeme jako „definici“ tuhosti.

J. D. Lambert definuje tuhost takto:

Pokud numerická metoda s konečnou oblastí absolutních stabilita, aplikován na systém s jakýmkoli počáteční podmínky, je nucen použít v určitém intervalu integrace délku kroku, která je příliš malá ve vztahu k plynulosti přesného řešení v tomto intervalu, pak se říká, že systém je tuhý v tom intervalu.

Existují i ​​další charakteristiky, které se projevují mnoha příklady tuhých problémů, ale pro každý z nich existují protiklady, takže tyto charakteristiky nedělají dobré definice tuhosti. Definice založené na těchto charakteristikách jsou nicméně některými autory běžně používány a jsou dobrým vodítkem ohledně přítomnosti tuhosti. Lambert je z výše uvedených důvodů označuje spíše jako „prohlášení“ než jako definice. Některé z nich jsou:

1. Systém lineárních konstantních koeficientů je tuhý, pokud je celý vlastní čísla mají zápornou skutečnou část a poměr tuhosti je velký.
2. K tuhosti dochází, když délka kroku omezují požadavky na stabilitu, nikoli požadavky na přesnost.
3. Tuhost nastává, když se některé složky roztoku rozpadají mnohem rychleji než jiné.^[3]

Etymologie

Původ pojmu „tuhost“ nebyl jednoznačně stanoven. Podle Joseph Oakland Hirschfelder se používá výraz „tuhý“, protože tyto systémy odpovídají těsnému spojení mezi řidičem a řízený v servomechanismy.^[4]Podle Richarda. L. Burden a J. Douglas Faires,

Při standardu mohou nastat značné potíže numerické techniky se používají k přiblížení řešení a diferenciální rovnice když přesné řešení obsahuje podmínky formuláře E^λt, kde λ je komplexní číslo se zápornou skutečnou částí.

...

Problémy s rychle se rozpadajícími přechodnými řešeními se přirozeně vyskytují v široké škále aplikací, včetně studia pružinových a tlumicích systémů, analýzy řídicí systémy a problémy v chemická kinetika. Toto jsou všechny příklady třídy problémů zvaných tuhé (matematické tuhosti) systémy diferenciálních rovnic, kvůli jejich aplikaci při analýze pohybu pružiny a hmoty systémy mít velké jarní konstanty (fyzický ztuhlost ).^[5]

Například problém počáteční hodnoty

${ displaystyle m { ddot {x}} + c { dot {x}} + kx = 0, qquad x (0) = x_ {0}, qquad { dot {x}} (0) = 0,}$

(10)

s m = 1, C = 1001, k = 1000, lze napsat ve tvaru (5) s n = 2 a

${ displaystyle mathbf {A} = left ({ begin {array} {rr} 0 & 1 - 1000 & -1001 end {array}} right),}$

(11)

${ displaystyle mathbf {f} (t) = left ({ begin {pole} {c} 0 0 end {pole}} right),}$

(12)

${ displaystyle mathbf {x} (0) = vlevo ({ začátek {pole} {c} x_ {0} 0 konec {pole}} doprava),}$

(13)

a má vlastní čísla${ displaystyle { overline { lambda}} = - 1 000, { podtržení { lambda}} = - 1}$. Obě vlastní čísla mají zápornou skutečnou část a poměr tuhosti je

${ displaystyle { frac {| -1000 |} {| -1 |}} = 1000,}$

(14)

což je docela velké. Systém (10) pak určitě splňuje výroky 1 a 3. Zde je jarní konstanta k je velké a tlumení konstantní C je ještě větší.^[6] (Všimněte si, že „velký“ je vágní, subjektivní výraz, ale čím větší jsou výše uvedené veličiny, tím výraznější bude účinek tuhosti.) Přesné řešení (10) je

${ displaystyle x (t) = x_ {0} left (- { frac {1} {999}} e ^ {- 1000t} + { frac {1000} {999}} e ^ {- t} vpravo) přibližně x_ {0} e ^ {- t}.}$

(15)

Všimněte si, že (15) se chová téměř jako jednoduchá exponenciál X₀E^−t, ale přítomnost E^−1000t i s malým koeficientem stačí k tomu, aby byl numerický výpočet velmi citlivý na velikost kroku. Stabilní integrace (10) vyžaduje velmi malou velikost kroku, dokud se nedostane do hladké části křivky řešení, což má za následek chybu mnohem menší, než je požadováno pro přesnost. Systém tedy splňuje i výrok 2 a Lambertovu definici.

A-stabilita

Chování numerických metod při tuhých problémech lze analyzovat použitím těchto metod na testovací rovnici y ' = ky podléhá původnímu stavu y(0) = 1 s ${ displaystyle k in mathbb {C}}$. Řešení této rovnice je y (t) = E^kt. Toto řešení se blíží nule jako ${ displaystyle t to infty}$ když ${ displaystyle mathrm {Re} , (k) <0.}$ Pokud numerická metoda také vykazuje toto chování (pro pevnou velikost kroku), pak se říká, že je A-stabilní.^[7] (Všimněte si, že numerická metoda, která je stabilní L (viz níže), má silnější vlastnost, že řešení se blíží nule v jediném kroku, protože velikost kroku jde do nekonečna.) A-stabilní metody nevykazují problémy s nestabilitou, jak je popsáno v motivační příklad.

Metody Runge – Kutta

Metody Runge – Kutta aplikován na testovací rovnici ${ displaystyle y '= k cdot y}$ mít formu ${ displaystyle y_ {n + 1} = phi (hk) cdot y_ {n}}$, a indukcí, ${ displaystyle y_ {n} = left ( phi (hk) right) ^ {n} cdot y_ {0}}$. Funkce ${ displaystyle phi}$ se nazývá funkce stability. Podmínka tedy, že ${ displaystyle y_ {n} až 0}$ tak jako ${ displaystyle n to infty}$ je ekvivalentní k ${ displaystyle | phi (hk) | <1}$. To motivuje k definici oblast absolutní stability (někdy označované jednoduše jako oblast stability), což je sada ${ displaystyle {z v mathbb {C} | , | phi (z) | <1 }}$. Metoda je A-stabilní, pokud oblast absolutní stability obsahuje sadu ${ displaystyle {z v mathbb {C} | mathrm {Re} (z) <0 }}$, tj. levá polovina roviny.

Příklad: Eulerovy metody

Růžový disk ukazuje oblast stability pro Eulerovu metodu.

Zvažte výše uvedené Eulerovy metody. Explicitní Eulerova metoda aplikován na testovací rovnici ${ displaystyle y '= k cdot y}$ je

${ displaystyle y_ {n + 1} = y_ {n} + h cdot f (t_ {n}, y_ {n}) = y_ {n} + h cdot (ky_ {n}) = y_ {n} + h cdot k cdot y_ {n} = (1 + h cdot k) y_ {n}.}$

Proto, ${ displaystyle y_ {n} = (1 + hk) ^ {n} cdot y_ {0}}$ s ${ displaystyle phi (z) = 1 + z}$. Oblast absolutní stability pro tuto metodu je tedy ${ displaystyle {z v mathbb {C} || 1 + z | <1 }}$ což je disk zobrazený vpravo. Eulerova metoda není stabilní.

Motivující příklad měl ${ displaystyle k = -15}$. Hodnota z při velikosti kroku ${ displaystyle h = 1/4}$ je ${ displaystyle z = -15 * 1/4 = -3,75}$, který je mimo oblast stability. Numerické výsledky ve skutečnosti nekonvergují k nule. S velikostí kroku ${ displaystyle h = 1/8}$, my máme ${ displaystyle z = -1,875}$ který je právě uvnitř oblasti stability a číselné výsledky konvergují k nule, i když poměrně pomalu.

Příklad: Lichoběžníková metoda

Růžová oblast je oblastí stability pro lichoběžníkovou metodu.

Zvažte lichoběžníkovou metodu

${ displaystyle y_ {n + 1} = y_ {n} + { tfrac {1} {2}} h cdot left (f (t_ {n}, y_ {n}) + f (t_ {n + 1}, y_ {n + 1}) right),}$

při aplikaci na testovací rovnici ${ displaystyle y '= k cdot y}$, je

${ displaystyle y_ {n + 1} = y_ {n} + { tfrac {1} {2}} h cdot left (ky_ {n} + ky_ {n + 1} right).}$

Řešení pro ${ displaystyle y_ {n + 1}}$ výnosy

${ displaystyle y_ {n + 1} = { frac {1 + { frac {1} {2}} hk} {1 - { frac {1} {2}} hk}} cdot y_ {n} .}$

Stabilitní funkce tedy je

${ displaystyle phi (z) = {1+ {1 nad 2} z nad 1- {1 nad 2} z}}$

a oblast absolutní stability je

${ displaystyle left {z in mathbb {C} left | left | {1+ {1 over 2} z over 1- {1 over 2} z} right | <1 dobře dobře}.}$

Tato oblast obsahuje levou polovinu roviny, takže lichoběžníková metoda je stabilní. Ve skutečnosti je oblast stability identická s levou polovinou roviny, a tedy numerické řešení ${ displaystyle y '= k cdot y}$ konverguje k nule, pokud a pouze pokud přesné řešení ano. Lichoběžníková metoda nicméně nemá dokonalé chování: tlumí všechny rozpadající se složky, ale rychle se rozpadající složky jsou tlumeny jen velmi mírně, protože ${ displaystyle phi (z) až 1}$ tak jako ${ displaystyle z to - infty}$. To vedlo ke konceptu L-stabilita: metoda je stabilní L, pokud je stabilní A a ${ displaystyle | phi (z) | až 0}$ tak jako ${ displaystyle z to infty}$. Lichoběžníková metoda je A-stabilní, ale ne L-stabilní. The implicitní Eulerova metoda je příkladem metody L-stabilní.^[8]

Obecná teorie

Stabilitní funkce a Metoda Runge – Kutta s koeficienty ${ displaystyle A}$ a ${ displaystyle b}$ darováno

${ displaystyle phi (z) = { frac { det (I-zA + zeb ^ {T})} { det (I-zA)}},}$

kde ${ displaystyle e}$ označuje vektor jednotkami. Tohle je racionální funkce (jeden polynomiální děleno jiným).

Explicitní metody Runge – Kutta mají a přísně nižší trojúhelníkový matice koeficientu ${ displaystyle A}$ a jejich stabilitní funkce je tedy polynom. Z toho vyplývá, že explicitní metody Runge – Kutta nemohou být stabilní.

Stabilitní funkce implicitních metod Runge – Kutta je často analyzována pomocí objednat hvězdy. Řádová hvězda pro metodu se stabilizační funkcí ${ displaystyle phi}$ je definována jako množina ${ displaystyle {z in mathbb {C} || phi (z) |> | mathrm {e} ^ {z} | }}$. Metoda je A-stabilní tehdy a jen tehdy, pokud její funkce stability nemá v levé rovině žádné póly a její řádová hvězda neobsahuje žádná čistě imaginární čísla.^[9]

Vícestupňové metody

Lineární vícestupňové metody mít formu

${ displaystyle y_ {n + 1} = součet _ {i = 0} ^ {s} a_ {i} y_ {ni} + h součet _ {j = -1} ^ {s} b_ {j} f (t_ {nj}, y_ {nj}).}$

Při aplikaci na testovací rovnici se stanou

${ displaystyle y_ {n + 1} = součet _ {i = 0} ^ {s} a_ {i} y_ {ni} + hk součet _ {j = -1} ^ {s} b_ {j} y_ {nj},}$

které lze zjednodušit na

${ displaystyle (1-b _ {- 1} z) y_ {n + 1} - součet _ {j = 0} ^ {s} (a_ {j} + b_ {j} z) y_ {nj} = 0 }$

kde z = hk. Toto je lineární relace opakování. Metoda je stabilní, pokud všechna řešení {y_n} relace opakování konverguje k nule, když Re z <0. Charakteristický polynom je

${ displaystyle Phi (z, w) = w ^ {s + 1} - součet _ {i = 0} ^ {s} a_ {i} w ^ {si} -z součet _ {j = -1 } ^ {s} b_ {j} w ^ {sj}.}$

Všechna řešení konvergují k nule pro danou hodnotu z pokud všechna řešení w z Φ (z,w) = 0 leží v jednotkovém kruhu.

Oblast absolutní stability pro vícestupňovou metodu výše uvedené formy je pak množinou všech ${ displaystyle z in mathbb {C}}$ pro které všechny w takové, že Φ (z,w) = 0 uspokojit |w| <1. Opět platí, že pokud tato sada obsahuje levou poloviční rovinu, je o vícestupňové metodě řečeno, že je A-stabilní.

Příklad: Metoda Adams – Bashforth druhého řádu

Růžová oblast je oblastí stability pro metodu Adams – Bashforth druhého řádu.

Pojďme určit oblast absolutní stability pro dvoustupňovou metodu Adams – Bashforth

${ displaystyle y_ {n + 1} = y_ {n} + h left ({ tfrac {3} {2}} f (t_ {n}, y_ {n}) - { tfrac {1} {2 }} f (t_ {n-1}, y_ {n-1}) vpravo).}$

Charakteristický polynom je

${ displaystyle Phi (w, z) = w ^ {2} - (1 + { tfrac {3} {2}} z) w + { tfrac {1} {2}} z = 0}$

který má kořeny

${ displaystyle w = { tfrac {1} {2}} { Big (} 1 + { tfrac {3} {2}} z pm { sqrt {1 + z + { tfrac {9} {4 }} z ^ {2}}} { Big)},}$

tedy oblast absolutní stability je

${ displaystyle left {z in mathbb {C} left | left | { tfrac {1} {2}} { Big (} 1 + { tfrac {3} {2}} z pm { sqrt {1 + z + { tfrac {9} {4}} z ^ {2}}} { Big)} doprava | <1 doprava. doprava }.}$

Tato oblast je zobrazena vpravo. Nezahrnuje celou levou polorovinu (ve skutečnosti zahrnuje pouze skutečnou osu mezi z = -1 a z = 0), takže metoda Adams – Bashforth není stabilní.

Obecná teorie

Explicitní vícestupňové metody nikdy nemohou být stabilní, stejně jako explicitní metody Runge – Kutta. Implicitní vícestupňové metody mohou být stabilní A, pouze pokud je jejich pořadí nejvýše 2. Druhý výsledek je známý jako druhý Dahlquist bariéra; omezuje užitečnost lineárních vícestupňových metod pro tuhé rovnice. Příkladem metody A druhého řádu druhého řádu je výše uvedené lichoběžníkové pravidlo, které lze také považovat za lineární vícestupňovou metodu.^[10]

Viz také

• Číslo podmínky
• Diferenciální inkluze, rozšíření pojmu diferenciální rovnice, které umožňuje diskontinuity, částečně jako způsob, jak se vyhnout některým problémům s tuhostí
• Explicitní a implicitní metody

Poznámky

Reference

• Burden, Richard L .; Faires, J. Douglas (1993), Numerická analýza (5. vydání), Boston: Prindle, Weber a Schmidt, ISBN  0-534-93219-3.
• Dahlquist, Germund (1963), "Speciální problém stability pro lineární vícestupňové metody", BIT, 3 (1): 27–43, doi:10.1007 / BF01963532, hdl:10338.dmlcz / 103497.
• Eberly, David (2008), Stabilitní analýza pro systémy diferenciálních rovnic (PDF).
• Ehle, B.L. (1969), Na Padé aproximace exponenciální funkce a A-stabilní metody pro numerické řešení problémů počáteční hodnoty (PDF), University of Waterloo.
• Gear, C. W. (1971), Numerické problémy počáteční hodnoty v obyčejných diferenciálních rovnicích, Englewoodské útesy: Prentice Hall.
• Gear, C. W. (1981), „Numerické řešení obyčejných diferenciálních rovnic: Je ještě co dělat?“, Recenze SIAM, 23 (1): 10–24, doi:10.1137/1023002.
• Hairer, Ernst; Wanner, Gerhard (1996), Řešení obyčejných diferenciálních rovnic II: Tuhé a diferenciálně-algebraické úlohy (druhé vydání), Berlín: Springer-Verlag, ISBN  978-3-540-60452-5.
• Hirshfelder, J. O. (1963), „Aplikovaná matematika používaná v teoretické chemii“, Symposium americké matematické společnosti: 367–376.
• Iserles, Arieh; Nørsett, Syvert (1991), Objednat hvězdy, Chapman & Hall, ISBN  978-0-412-35260-7.
• Kreyszig, Erwin (1972), Pokročilá inženýrská matematika (3. vyd.), New York: Wiley, ISBN  0-471-50728-8.
• Lambert, J. D. (1977), D. Jacobs (ed.), „Problém počáteční hodnoty pro obyčejné diferenciální rovnice“, Současný stav v numerické analýze, New York: Akademický tisk: 451–501.
• Lambert, J. D. (1992), Numerické metody pro běžné diferenciální systémy, New York: Wiley, ISBN  978-0-471-92990-1.
• Mathews, John; Fink, Kurtis (1992), Numerické metody využívající MATLAB.
• Stiskněte, WH; Teukolsky, SA; Vetterling, WT; Flannery, BP (2007). „Oddíl 17.5. Tuhé sady rovnic“. Numerické recepty: Umění vědecké práce na počítači (3. vyd.). New York: Cambridge University Press. ISBN  978-0-521-88068-8.
• Shampine, L. F .; Gear, C. W. (1979), „Pohled uživatele na řešení tuhých obyčejných diferenciálních rovnic“, Recenze SIAM, 21 (1): 1–17, doi:10.1137/1021001.
• Wanner, Gerhard; Hairer, Ernst; Nørsett, Syvert (1978), „Pořadí hvězd a teorie stability“, BIT, 18 (4): 475–489, doi:10.1007 / BF01932026.
• Stabilita metod Runge-Kutta [1]

externí odkazy

• Úvod do fyzikálně založeného modelování: energetické funkce a tuhost
• Tuhé systémy Lawrence F. Shampine a Přeskočit Thompsona Scholarpedia, 2(3):2855. doi: 10,4249 / scholarpedia.2855