Paprastųjų diferencialinių lygčių skaitiniai sprendiniai

Apibrėžkite modelį ir laiko tinklelį

Pateikto padėties modelio pradinė vertė yra 1 m, o kitimo greitis intervale [0 s, 2 s] yra x′(t) = (2 m/s²)t. Ar galime apskaičiuoti jo būseną trumpais laiko žingsniais? Čia t yra praėjęs laikas, o x – padėtis su ženklu, išreikšta metrais išilgai fiksuotos ašies, stebint tą patį tašką ir nekeičiant atskaitos pradžios. Tikslusis sprendinys, naudojamas patikrai, yra x(t) = 1 m + (1 m/s²)t². Kai laikas yra 2 s, jis duoda 5 m; palyginsime su šia verte skaitinius rezultatus, gautus naudojant du skirtingus žingsnio dydžius.

Brūkšnelis žymi išvestinę laiko atžvilgiu. Bendresniu atveju skaliariniame pradinės vertės uždavinyje nurodomas kitimo greitis y′(t) = F(t, y(t)) ir pradinė vertė y(t₀) = y₀. Reikia suprasti išvestinės ir pradinės sąlygos reikšmę; tam parengia paskelbtos pamokos „Išvestinės, Teiloro eilutės ir vietinis artinys“ bei „Paprastosios diferencialinės lygtys“. Čia naudojamas padėties modelis yra sutartinis bandymo pavyzdys su žinomu sprendiniu, o ne išmatuota trajektorija ar patikrintas jėgos dėsnis. Kituose uždaviniuose darykite prielaidą, kad sprendinys egzistuoja numatytame intervale; skaitinis metodas šio egzistavimo neįrodo.

Pasirinkite teigiamą sveikąjį žingsnių skaičių N, žingsnio dydį h > 0 ir laiko momentus tₙ = t₀ + nh, kai n = 0, …, N. Apatinis indeksas n skaičiuoja žingsnius. Nustatykite h = (t_N − t₀)/N, kad paskutinis tinklelio taškas atitiktų pasirinktą galutinį laiko momentą. Apskaičiuota vertė Yₙ yra y(tₙ) artinys; ji nėra nei naujas matavimo prietaiso rodmuo, nei visas tolydus sprendinys tarp tinklelio taškų.

Žengimas išreikštiniu Eulerio metodu

Išreikštiniame Eulerio metode naudojama formulė Yₙ₊₁ = Yₙ + hF(tₙ, Yₙ), pradedant nuo Y₀ = y₀. Atliekant šį atnaujinimą, kitimo sparta žingsnio pradžioje laikoma pastovia visą žingsnį. Jei y yra padėtis, F matuojama m/s, o hF – m, todėl sumos dydžio rūšis ir matavimo vienetai yra tinkami.

Kitimo spartą reikia apskaičiuoti pagal dabartinę skaitinę būseną, o ne pagal nežinomą tikslią būsimą būseną. Mūsų padėties modelyje F nepriklauso nuo padėties, todėl Xₙ₊₁ = Xₙ + h(2 m/s²)tₙ. Kai h = 1 s, iš X₀ = 1 m gauname X₁ = 1 m ir X₂ = 3 m. Tikslios padėtys tais laiko momentais yra 2 m ir 5 m.

Net ir tiksliai laikantis atnaujinimo taisyklės, gauti skaičiai nėra tikslūs diferencialinės lygties sprendiniai. Po pirmojo Eulerio metodo žingsnio padėtis lieka nepakitusi, nes pradinė kitimo sparta lygi nuliui, nors per tą žingsnį nurodyta kitimo sparta tampa teigiama. Kaip interpoliuoti vaizdą tarp taškų, pasirenkama atskirai; šioje pamokoje lyginamos reikšmės tais pačiais laiko tinklelio momentais.

Vieno žingsnio ir globaliosios paklaidos atskyrimas

Norėdami atskirai įvertinti vieno žingsnio aproksimacijos paklaidą, įsivaizduokime, kad tą žingsnį pradedame nuo tikslios y(tₙ) reikšmės ir skaičiuojame tiksliai. Apibrėžkime dₙ = y(tₙ₊₁) − [y(tₙ) + hF(tₙ, y(tₙ))]. Tai vieno žingsnio būsenos neatitiktis, kurios matavimo vienetai tokie pat kaip y. Terminija įvairuoja: kai kuriuose šaltiniuose nukirtimo arba suderinamumo paklaida vadinamas dydis dₙ/h, todėl svarbu atsižvelgti į formulę ir matavimo vienetus.

Jei tiksliojo sprendinio antroji išvestinė yra tolydi visame žingsnio intervale, pagal Teiloro teoremą dₙ = (h²/2)y″(ξ), kur ξ yra tam tikras vidinis šio intervalo laiko taškas. Todėl visame intervale galiojantis įvertis |y″| ≤ M duoda |dₙ| ≤ Mh²/2. Jei y yra padėtis, M matavimo vienetai yra m/s². Šis įvertis gaunamas pradedant nuo tikslios būsenos; jis neapima paklaidos, paveldėtos iš ankstesnių skaitinio skaičiavimo žingsnių.

Globalioji paklaida tinklelio taške yra eₙ = y(tₙ) − Yₙ, o jos absoliutusis dydis yra Eₙ = |eₙ|. Praktiškai Eulerio metodo žingsniai paprastai pradedami nuo apytikslės Yₙ reikšmės. Pasikeitusi būsena gali pakeisti kito žingsnio kitimo spartą, todėl globaliajai paklaidai įvertinti reikia atsižvelgti ir į naujas žingsnių neatitiktis, ir į paveldėtas paklaidas. Apvalinimas, netikslūs pradiniai duomenys ir klaidingas fizikinis modelis yra papildomi paklaidų šaltiniai.

Paklaidos sklidimas esant aiškiai nurodytoms sąlygoms

Tarkime, kad kiekvienu tinklelio laiko momentu lyginamoms būsenoms galioja |F(t,u) − F(t,v)| ≤ L|u − v|, kur L ≥ 0. Tai Lipšico įvertis: jis apriboja, kiek gali pasikeisti kitimo sparta pasikeitus būsenai. L matavimo vienetai yra s⁻¹. Reikalaukime, kad tikslusis sprendinys ir apskaičiuotos būsenos liktų srityje, kurioje galioja šis įvertis, o antroji išvestinė būtų tolydi ir tenkintų |y″| ≤ M visame skaičiavimo intervale.

Iš tiksliojo žingsnio tapatybės atėmę Eulerio metodo atnaujinimo formulę ir paėmę absoliučiąsias vertes, gauname Eₙ₊₁ ≤ (1 + hL)Eₙ + Mh²/2. Daugiklis 1 + hL yra bematis; kiekvieno dėmens matavimo vienetai tokie pat kaip y. Tai viršutinis įvertis, o ne teiginys, kad paklaidos visada yra vienodo ženklo ar visada didėja. Jei neįtraukiami papildomi paklaidos dėmenys, daroma prielaida, kad skaičiuojama tiksliai ir naudojami nurodyti pradiniai duomenys.

Mūsų padėties kitimo greičiui L = 0 s⁻¹ ir M = 2 m/s². Kai X₀ tikslus, pakartotinai taikant įvertį gaunama E_N ≤ N Mh²/2. Fiksuotu galutiniu laiko momentu 2 s tai yra (1 m/s²)(2 s)h: 2 m, kai h = 1 s, ir 1 m, kai h = 0.5 s. Ši konkreti kvadratinė funkcija pasiekia šias įverčių ribas. Bendruoju atveju nuo būsenos priklausantiems kitimo greičiams reikia atskiro paklaidos sklidimo įverčio; vien mažesnis lokalusis neatitikimas nėra išsami globaliosios paklaidos analizė.

Įvertinkite du kitimo greičius išreikštiniame žingsnyje

Taikant dviejų stadijų alternatyvą, pirmiausia apskaičiuojama k₁ = F(tₙ, Yₙ) ir prognozuojama reikšmė P = Yₙ + hk₁. Tada apskaičiuojama k₂ = F(tₙ + h, P) ir priimama Yₙ₊₁ = Yₙ + (h/2)(k₁ + k₂). Abi k reikšmės yra kitimo greičiai. Ši prognozavimo ir vidurkinimo formulė yra išreikštinis patobulinto Eulerio metodo variantas; P yra bandomoji būsena, o ne priimta kito žingsnio reikšmė.

Mūsų padėties modelyje, kai t₀ = 0 s ir h = 1 s, k₁ = 0 m/s, P = 1 m ir k₂ = 2 m/s. Priimta X₁ reikšmė yra 2 m. Kartojant nuo 1 s, apskaičiuojamas 2 m/s ir 4 m/s vidurkis ir gaunama X₂ = 5 m. Kitimo greičiai yra afininės laiko funkcijos ir nepriklauso nuo būsenos, todėl jų reikšmių žingsnio galuose vidurkis duoda tikslų padėties pokytį kiekviename šio pateikto modelio žingsnyje.

Šis tikslumas būdingas tik šiai kitimo greičio taisyklei. Atskirai nagrinėjant mažėjimo greitį F(t,y) = −λy, tos pačios dviejų stadijų formulės daugiklis yra 1 − λh + (λh)²/2. Kai λh = 2, jis lygus 1, todėl atnaujinta reikšmė nemažėja, nors kitimo greitis įvertinamas du kartus. Papildomos stadijos, metodo pavadinimas ar tikslioji aritmetika savaime neužtikrina reikiamo tikslumo ar stabilumo.

Patikrinkite mažėjimo stabilumą ir ženklus

Naudokite atskirą apibrėžtą modelį: θ(t) yra temperatūros perviršis kelvinais, palyginti su fiksuota 280 K atskaitos temperatūra. Duota θ′ = −λθ, λ = 2 s⁻¹ ir θ(0 s) = 8 K intervale [0 s, 3 s]. Tikslusis sprendinys θ(t) = (8 K)e^(−λt) šiame intervale yra teigiamas ir mažėja; laipsnio rodiklis yra bedimensis. Tai mokomasis kitimo greičio dėsnis, o ne realaus vėstančio kūno modelio patikra.

Išreikštinis Eulerio metodas duoda Θₙ₊₁ = rΘₙ, kur r = 1 − λh. Pradinės reikšmės pokytis kiekviename žingsnyje dauginamas iš r. Šioje tiesinėje bandomojoje rekurentinėje formulėje |r| < 1 slopina šio pokyčio modulį, o |r| > 1 jį didina. Griežto slopinimo sąlyga yra 0 < λh < 2. Kai h = 1 s, r = −1: 8 K, −8 K, 8 K rodo ženklų kaitą be slopinimo. Kai h = 1.5 s, r = −2 duoda 8 K, −16 K, 32 K iki 3 s.

Slopinimas nėra tas pats, kas ženklo išsaugojimas ir reikšmių tikslumas. Kai h = 0.75 s, r = −0.5 slopina modulį, tačiau pirmoji reikšmė yra −4 K, priešingai nei tikslusis teigiamas perviršis. Kad Eulerio metodo reikšmės kiekviename žingsnyje, pradedant nuo teigiamos pradinės reikšmės, būtų neneigiamos, turi galioti 0 ≤ λh ≤ 1; kai λh = 1, jos iš karto tampa lygios nuliui, užuot sutapusios su teigiamu tiksliuoju sprendiniu. Šie bandymai taikomi šiai skaliarinei rekurentinei formulei ir nėra universalus netiesinių ar susietųjų sistemų kriterijus.

Atpažinkite neišreikštinį atnaujinimą

Neišreikštinis Eulerio metodas apibrėžiamas lygtimi Yₙ₊₁ = Yₙ + hF(tₙ₊₁, Yₙ₊₁). Nežinoma kito žingsnio būsena įeina į jos pačios kitimo greičio skaičiavimą, todėl žingsniui atlikti paprastai reikia išspręsti lygtį. Netiesiniam žingsniui gali reikėti papildomų prielaidų ir lygčių sprendimo algoritmo; pavadinimas „neišreikštinis“ negarantuoja, kad ši lygtis turi vienintelį sprendinį.

Mūsų tiesinio mažėjimo uždavinyje pertvarkius gaunama Θₙ₊₁ = Θₙ/(1 + λh). Kadangi λ > 0 ir h > 0, vardiklis yra teigiamas, o daugiklis yra griežtai tarp nulio ir vieneto. Šis konkretus atnaujinimas išsaugo teigiamą perviršį ir slopina pradinės reikšmės trikdžius esant bet kokiam teigiamam žingsniui, kuris neperžengia nurodyto skaičiavimo intervalo.

Kai h = 1 s, pirmoji neišreikštiniu metodu gauta vertė yra 8/3 K; tiksli vertė yra (8 K)e^(−2), ir ji skiriasi. Taigi skaitinis slopinimas negarantuoja tikslumo. Šis palankus daugiklis yra pateikto tiesinio bandomojo uždavinio savybė, o ne teiginys, kad kiekvienas neišreikštinis metodas yra stabilus ar tiksliai atspindi fiziką sprendžiant bet kurią lygtį.

Patikrinkite skaičiavimą ir modelį

Užrašykite lygtį, pradinius duomenis, matavimo vienetus, būsenų sritį, laiko intervalą, metodą, žingsnių dydžius ir aritmetinių skaičiavimų tikslumą. Nepriklausomai patikrinkite kiekvieną kitimo spartos skaičiavimą ir vertės atnaujinimą, o skirtingais žingsnių dydžiais gautus rezultatus palyginkite tuo pačiu galutiniu laiko momentu. Mūsų padėties modeliui Eulerio metodu laiko momentu 2 s gaunamos vertės 3 m, 4 m ir 4.5 m, kai žingsnių dydžiai yra atitinkamai 1 s, 0.5 s ir 0.25 s; tiksli kontrolinė vertė yra 5 m.

Žingsnio smulkinimas yra naudingas patikrinimas, tačiau vien dviejų skaičiavimų rezultatų sutapimas nesuteikia paklaidos garantijos. Paklaidos rėžiui būtinos atitinkamos glodumo, jautrumo būsenai ir srities sąlygos. Daugiau žingsnių taip pat reiškia daugiau aritmetinių operacijų, todėl mažinant h riboto tikslumo skaičiavimuose bendroji paklaida nebūtinai mažėja. Nulinis neatitikimas atnaujinant skaitinę vertę tik rodo, kad buvo laikomasi pasirinktos rekurentinės formulės.

Tikslus pateiktos lygties sprendinys tebėra to modelio rezultatas. Tikrinant realią sistemą taip pat reikia duomenų, pagrindžiančių jos kitimo spartos dėsnį, parametrus, pradinius duomenis ir veikimo sąlygas. Laikykitės nurodytų intervalų net tada, kai algebrinių formulių vertes galima apskaičiuoti ir už jų ribų. Kitoje pamokoje „Skaitiniai dalinių išvestinių diferencialinių lygčių sprendiniai“ skaitinis modeliavimas papildomas erdviniu kitimu ir kraštinėmis sąlygomis.

Išnagrinėti pavyzdžiai

Pabandykite patys

1. Padėties modeliui atlikite du išreikštinio Eulerio metodo žingsnius, kai h = 0.5 s, nuo 0 s iki 1 s. Raskite X₁ ir X₂, tada palyginkite X₂ su x(1 s).

   Rodyti atsakymą ir paaiškinimą

   X₁ = 1 m ir X₂ = 1.5 m. Tiksli padėtis, su kuria lyginama, yra 2 m, todėl e₂ = +0.5 m.

   Greičiai žingsnių pradžioje yra 0 m/s ir 1 m/s, todėl prieaugiai yra 0 m ir 0.5 m. Tai, kad pirmojo žingsnio pradžioje greitis lygus nuliui, nereiškia, kad tolydusis greitis lygus nuliui visame tame intervale. Abu laiko taškai priklauso nurodytai apibrėžimo sričiai.

2. Tarkime, padėtį nusakantis sprendinys turi tolydžią antrąją išvestinę, o 0.2 s trukmės žingsnyje galioja |y″| ≤ 6 m/s². Koks |dₙ| įvertis galioja, jei pradedama nuo tikslios būsenos ir skaičiuojama tiksliai? Ar jis savaime yra per ilgą skaičiavimą sukauptos paklaidos įvertis?

   Rodyti atsakymą ir paaiškinimą

   Vieno žingsnio būsenos paklaidos viršutinis įvertis yra 0.12 m. Jis savaime nėra globaliosios paklaidos įvertis.

   Mh²/2 = (6 m/s²)(0.2 s)²/2 = 0.12 m. Pagal prielaidą žingsnis iš naujo pradedamas nuo tikslios būsenos. Kai žingsnių seka pradedama nuo ankstesnių apytikslių reikšmių, taip pat reikia atsižvelgti į paveldėtą paklaidą ir nuo būsenos priklausančius greičio pokyčius.

3. Taikykite pateiktas Eulerio metodo paklaidos rekurentinio sąryšio sąlygas, kai L = 3 s⁻¹, h = 0.2 s, M = 6 m/s² ir Eₙ ≤ 0.1 m. Koks iš jų gaunamas Eₙ₊₁ viršutinis įvertis?

   Rodyti atsakymą ir paaiškinimą

   Eₙ₊₁ ≤ 0.28 m, jei tenkinamos pateiktos reguliarumo, būsenų srities ir tikslios aritmetikos sąlygos.

   1 + hL = 1.6. Paveldėtos paklaidos narys neviršija 1.6(0.1 m) = 0.16 m, o naujasis vieno žingsnio narys neviršija 0.12 m. Jų suma yra 0.28 m. Šis viršutinis įvertis nereiškia, kad tikroji paklaida turi būti jam lygi ar viršyti ankstesnę paklaidą.

4. Kai λ = 2 s⁻¹, o pradinis perviršis yra 8 K, išreikštiniame Eulerio metode naudojamas h = 0.75 s žingsnis. Raskite metodo daugiklį ir pirmąją reikšmę. Ar slopimas garantuoja, kad išliks teigiamas tikslaus sprendinio ženklas?

   Rodyti atsakymą ir paaiškinimą

   r = −0.5, o Θ₁ = −4 K. Jos modulis mažėja, tačiau ženklas neišlieka.

   λh = 1.5 patenka į griežto slopinimo intervalą (0,2). Neigiamas daugiklis pakeičia ženklą. Tikslus temperatūros perviršis (8 K)e^(−1.5) yra teigiamas. Čia lyginamas temperatūros perviršis, o ne teigiama, kad absoliučioji temperatūra yra neigiama.

5. Mažėjimo modeliui taikykite dviejų etapų prognozavimo ir vidurkinimo žingsnį, kai h = 1 s. Raskite prognozuojamąją vertę, priimtą vertę ir daugiklį. Ar gaunamas griežtas mažėjimas?

   Rodyti atsakymą ir paaiškinimą

   P = −8 K, priimtas temperatūros perviršis yra 8 K, o jo daugiklis yra 1. Griežtas mažėjimas negaunamas.

   Kai k₁ = −16 K/s, gauname P = 8 K + (1 s)(−16 K/s) = −8 K. Tuomet k₂ = +16 K/s. Skaičiuojant priimtą vertę pridedama pusė sekundės, padauginta iš nuliui lygios šių greičių sumos, todėl lieka 8 K. Daugiklio formulė, kai λh = 2, taip pat duoda 1. Du kitimo greičio įvertinimai šio stabilumo apribojimo nepašalina.

6. Kai h = 1 s, padėties vertės 1 m, 1 m, 3 m tiksliai tenkina Eulerio rekurentinę formulę. Ar tai įrodo, kad jos yra tikslūs pradinės paprastosios diferencialinės lygties sprendiniai arba kad realus objektas juda pagal šią lygtį? Kokia patikra čia gali padėti atskirti skaitines ir tikslias vertes?

   Rodyti atsakymą ir paaiškinimą

   Nė viena iš šių išvadų neišplaukia. Įstačius į žinomą patikrai skirtą sprendinį, vėlesniuose tinklelio taškuose gaunamos 1 m ir 2 m paklaidos.

   Nulinis žingsnio skaičiavimo neatitikimas patvirtina, kad laikomasi nurodytos skaitinės rekurentinės formulės. Nurodytas tolydusis sprendinys laiko momentais 1 s ir 2 s duoda 2 m ir 5 m, o ne apskaičiuotas 1 m ir 3 m vertes. Smulkinant tinklelį ir išlaikant tą patį galutinį laiko momentą, galima atlikti dar vieną skaitinę patikrą; fiziniam pagrįstumui nustatyti reikėtų nepriklausomų stebėjimų ir modelį pagrindžiančių duomenų.

Ką verta įsidėmėti

• Paprastųjų diferencialinių lygčių sprendimo algoritmas apytiksliai apskaičiuoja pradinės vertės uždavinio sprendinį nurodytais laiko momentais; jo gautos vertės nėra nauji matavimai.
• Išreikštinis Eulerio metodas naudoja esamą kitimo greitį, o kiekviename skaičiavimo žingsnyje dydžių vienetai turi išlikti suderinti.
• Vieno žingsnio būsenos paklaida nustatoma darant prielaidą, kad pradinė būsena yra tiksli; globalioji paklaida taip pat apima ankstesnių paklaidų poveikį.
• Paklaidų sklidimui nagrinėti reikia glodumo, jautrumo būsenai ir apibrėžimo srities prielaidų, o kitiems paklaidų šaltiniams įvertinti – papildomų narių.
• Slopinimas, reikšmių teigiamumas ir jų tikslumas yra skirtingos savybės net ir paprastame tiesiniame mažėjimo modelyje.
• Prieš interpretuodami apskaičiuotą trajektoriją, patikrinkite žingsnio smulkinimo rezultatus tuo pačiu laiko momentu, tikslius etaloninius atvejus ir modelį pagrindžiančius duomenis.