Перейти к содержанию

Численные методы решения задачи Коши для обыкновенных дифференциальных уравнений

Материал из Викиверситета

Найти функцию для y(t) удовлетворяющую уравнению

Задача Коши для ОДУ
{y′=f(t,y(t)),t>t0y(t0)=y0

Ищем интегральную кривую, проходящую через (t0,y0)

Применяя формулу численного дифференцирования

[править]
y(ti+1)−y(ti)h=f(ti,y(ti))⇒yi+1=yi+h⋅f(ti,yi)

Формула Тейлора

[править]

Пусть y(t) - найдено. Найдем t+h:

y(t+h)=y(t)+hy′(t)+h22y″(t)+...+hpp!y(p)(t)+hp+1y(p+1)(ξ)(p+1)!

Первая производная - это правая часть.

II-я: y″=f't+f'yy't=f't+f'yf

III-я: y‴=f′'tt+f′'tyf+(f'yt+f'yy)+f'y(f't+f'yf)

  • Формула Тейлора I-го порядка точности по h (Метод Эйлера)
y(t+h)=y(t)+hf(t,y(t))
  • Формула Тейлора II-го порядка точности по h:
y(t+h)=y(t)+hf(t,y(t))+12h2(f't+f'yf)

Метод Эйлера

[править]
{y′=f(t,y(t)),ti=t0+ih,i=1,n‾y(t0)=y0

Дискретизация задачи (переход от непрерывной функции к дискретной)

y(ti) - точное решение. yi - приближенное. y′(ti)≈h−1(yi+1−yi), yi+1−yih=f(ti,yi)

yn=(y0⋮yN) - в ответе получаем вектор.

{yi+1=yi+hf(ti,yi)y0 - задано шаг обычно берут h = 0.1}

Геометрическая интерпретация.

Геометрическая интерпретация метода Эйлера

l(t)=y(ti)+y′(ti)(t−ti) - уравнение касательной.

Погрешности:

Определение. Погрешность на шаге определяется по формуле ri=y(ti+1)−yi+1, где y(ti) - точное решение, а yi - приближенное.


Определение. Локальная погрешность: li=y(ti)−yi, если yi−1≡y(ti−1) - точное.


Определение. Глобальная погрешность: E(h)=max⁡0≤i≤N−1‖ri‖


Определение. Будем называть метод p-го порядка точности по h, если E(h)≤chp, где c - константа.


Модификации метода Эйлера

[править]

yi+1=yi+hki. изменяя ki можно легко получить меньшую погрешность.

  • ki=f(ti+12,yi+12) - усовершенствованный метод Эйлера.
{yi+12=yi+h2f(ti,yi)yi+1=yi+hf(ti+12,yi+12)⇔yi+1=yi+hki,ki=f(ti+12,yi+h2f(ti,yi)) (за угол наклона берется наклон касательной в средней точке).
E(h)=ch2
  • Метод Эйлера-Коши:

yi+1+hki,ki=ki1+ki22,ki1=f(ti,yi),ki2=f(ti+h,yi+hki2). E(h)=ch2.

Методы Рунге-Кутты

[править]

Метод Эйлера - это метод Рунге-Кутты I-го порядка точности. yi+1=yi+hki,ki - угловой коэффициент секущей.

m-этапный метод Рунге-Кутты

yi=y(ti)

Метод Рунге-Кутты

ti(1)=ti ti(2)=ti(1)+α2h0<α2<1ti(m)=ti(1)+αmh0<αm<1

ki(1)=f(ti(1),yi(1))=f(ti,yi)

ki(2)=f(ti(1)+α2h,yi+hβ21ki(1))

ki(3)=f(ti(1)+α3h,yi+h(β31ki(1)+β32ki(2)))

…

ki(m)=f(ti(1)+αmh,yi+h(βm1ki(1)+...+βm,m−1ki(m−1)))

Все коэффициенты (α и β) подлежат определению: ki=c1ki(1)+c2ki(2)+...+cmki(m) коэффициенты ci также неизвестны.

y(t+h)=y(t)+hy′(t)+h22y″(t)+...+hpp!y(p)(t)+O‾‾(hp+1)

Если хотим построить метод p-го порядка точности по h, yi+1=yi+hki. Возьмем разложение для точного решения. Правую часть раскладываем в ti. Коэффициенты α,β и c подбираем так, чтобы разложение совпало до p-го порядка включительно. Построим параметрическое семейство двухэтапных методов Рунге-Кутты (2-го порядка точности по h)

ti(1)=ti,ti(2)=ii+αh

ki(1)=f(ti,yi),ki(2)=f(ti+αh,yi+hβki(1),ki=c1ki(1)+c2ki(2)

yi+1=yi+hki=yi+h[c1f(ti,yi)+c2f(ti+αh,yi+hβki(1))]

f(ti+αh,yi+hβki(1))=f(ti,yi)+f't(ti,yi)αhf'y(ti,yi)hβki(1)+O‾‾(h2)

yi+1=yi+h[c1f(ti,yi)+c2(f(ti,yi)+f't(ti,yi)αh+hβf'y(ti,yi)ki(1))]+O‾‾(h3) - разложение приблзительного решения.

В полученное разложение будем подставлять точное решение, будем считать α=β

y″(ti)=f't(ti,y(ti))+f'y(t−i,y(ti))f(ti,yi(ti))

y(ti+1)=y(ti)+h(c1+c2)f(ti,y(ti))+c2αh2[f't(ti,y(ti))+f'y(ti,y(ti))f(ti,y(ti))]⇒y(ti+1)=y(ti)+h(c1+c2)y′+αc2h2y″(ti)+O‾‾(h3)

{c1+c2=1c2α=12⇒{c2=12αc1=1−12α,α - параметр

Имеем:

{yi+1=y+hkiki(1)=f(ti,yi)ki(2)=f(ti+αh,yi+αhki(1))ki=(1−12α)ki(1)+12αki(2)

При α=1 получаем метод Эйлера-Коши, α=12 - усовершенствованный метод Эйлера.

По построению: локальная погрешность этих методов для p-го порядка точности l=chp+1

Каноническая форма записи методов численного интегрирования

[править]

yi+1=yi+hf(ti,yi)→yi+1−yih=f(ti,yi) - явный одношаговый метод. 1h∑j=0kαjyi+1−j=Φ(ti,yi+1−k,...,yi,yi+1,h),α0≠0 - канонический вид формулы численного интегрирования. (неявный k-шаговый метод)

Если правая часть этой формулы не содержит yi+1, тогда это явный k-шаговый метод. Все методы Рунге-Кутты - явные, одношаговые.

yi+1−yif=Φ(ti,yi,h) - общая формула методов Р-К (каноническая формула явных одношаговых методов).

Определение. Аппроксимацией будем называть сеточную функцию Ψih, определенную так: \\Ψih=1h∑j=0kαiy(ti+1−j)−Φ(ti,y(ti+1−k),...,y(ti),y(ti+1),h)


Определение. Будем говорить, что метод сходится с порядком p, если величина E(h)=max⁡0≤1≤N|yi−y(ti)|≤chp


Погрешность аппроксимации метода Эйлера: Ψih=h22y″(t)

Устойчивость ЗК

[править]

{y′=f(t,y)y(t0)=y0{y~′=f(t,y~)−ϕy~(t0)=y~0 - возмущенная задача

ε(t)=y(t)−y~(t)⇒{ε′=f(t,y)−f(t,y~)+ϕε(t0)=y0−y~0=ε0

f(t,y)−f(t,y~)=f'y(t,y~~)(y−y~)⇒{ε′=λε+ϕ,λ=f'y(t,y~~),y~~=[y,y~]ε(t0)=ε0

Получим линейное неоднородное уравнение

{ε'1=ϕε(t0)=0 и {ε'2=λε2ε2(t0)=ε0⇒ε1=∫t0tϕ(τ)dτ,ε2=ε0exp⁡{∫t0tλ(τ)dτ}

Если λ(t)>0⇒ первое слагаемое ε2→[t→∞]∞. Если λ(t)<0⇒ε2→[t→∞]0

Оценка устойчивости ЗК

|ε(t)|≤k(τ)(|ε0|+∫t0τ|ϕ(t)|dt),k(τ)={exp⁡{∫t0τλ(t)dt},λ(t)≥01,λ(t)<0

Из этого неравенства вытекает устойчивость ЗК.

Рассмотрим возмущенную дискретную ЗК:

1h∑j=0ky~i+1−jλj=Φ(ti,y~i+1−k,...,y~i,y~i+1,h)−ϕ

ε0,ε1,...,εk−1 - погрешности вычисления стартовых точек метода.

Определение. Будем называть метод (1) устойчивым, если:

max⁡0≤i≤n|yi−y~i|≤k(τ)max⁡0≤i≤k−1(|εi|+∑i+kN−1h|ϕ|)


Многошаговые методы. Методы Адамса

[править]

Приближение правой части и интегрирование в многошаговых методах

∫titi+1y′dt=∫titi+1f(t,y(t))dt. y(t) приближаем полиномом по известным значениям.

Метод Адамса

Формулы Адамса-Башфорта

Метод Адамса-Башфорта

p1(t)=fi−1+fi−fi−1h(t−ti−1) - интерполяционный многочлен Ньютона. Интегрируемые функции f(t,y): ∫titi+1p1(t)dt=∫titi+1(fi−1+fi−fi−1h(t−ti−1))dt=fi−1h+fi−fi−1h(t−ti−1)22|t1ti+1=h2(3fi−fi−1)⇒yi+1=yi+h2(2fi−fi−1) - формула А-Б (двухшаговый метод, II порядка точности по h) Формулы Адамаса-Мултона

Метод Адамса-Моултона

Q1(t)=fi+fi+1−fih(t−ti)

∫titi+1Q1(t)dt=fih+(fi+1−fi)(t−ti)22h|t1ti+1=h2(fi+fi+1)

yi+1=yi+h2(fi+fi+1) - II порядка точности, более устойчив.

Порядок аппроксимации в построенных методах

[править]

Ψi - дискретная функция. Канонический вид: yi+1−yih=12fi+12fi+1

y′=f(t,y(t))fi=y′(ti)fi+1=y′(ti+1)

Ψi=y(ti+1)−y(ti)h−12y′(ti)−12y′(ti+1)=y(ti)+hy′(ti)h22y″(ti)+h36y‴(ti)+O__(h4)−y(t1)h−12y′(ti)−12(y′(ti)+hy″(ti)+h22y‴(ti)+O__(h3))=h2y″(ti)+h6y‴(ti)+O__(h3)−h2y″(ti)−h24y‴(ti)+O__(h3)=−h212y‴(ti)+O__(h3)

Ψi=−h212y‴(ti),Ψ=max⁡0≤i≤N|Ψi|≤M3h212⇒ II аппроксимационный порядок

Если функция y(t) - аналитична и k раз дифференцируема, то k - шаговый метод А-Б и k−1 шаговый метод А-М имеют k-й порядок аппроксимации по h.

Каноническая формула метода:

(A):1h∑j=0kαjyi+1−j=Φ(ti,yi+1−k,...,yi,yi+1) - дискретная ЗК. Стартовые точки y0,y1,...,yk−1

(B):1h∑j=0kαjyi+1−j∗=Φ(ti,yi+1−k∗,...,yi∗,yi+1∗) - дискретная ЗК с возмущенными данными. Стартовые точки y0,y1,...,yk−1

Определение. (Устойчивость k-шагового метода)
max⁡k≤i≤N|yi−yi∗|≤c{max⁡0≤i≤k−1|yi−yi∗|+max⁡|∑i=kNhΨi|}


Теорема. Если численный метод устойчив и имеет порядок аппроксимации p, то он сходится с p-м порядком точности по h.
Доказательство. Подставим в численный метод точное решение (A):

y(ti):(A)→(B),Ψi - погрешность аппроксимации Ψi=chp

max⁡k≤i≤N|yi−yih|≤c{max⁡0≤i≤k−1|yi−y(ti)|+max⁡|∑i=kNhc1hp|}=c{max⁡|yi−y(ti)|+c(T−t0)hp}

Предположим: max⁡k≤i≤N|yi−yih|≤c3hp, где c3=c(c2+c1(T−t0))

Если |yi−y(ti)|≤c(i)h0, то c2=max⁡c(i)


Уравнения высокого порядка и системы ДУ

[править]

{y'1=f1(t,y,...,ym)y'2=f2(t,y,...,ym)⋮y'm=fm(t,y,...,ym)y'1(t0)=y10⋮y'm(t0)=ym0 - нелинейная система ОДУ 1-го

y=(y1⋮ym)y0=(y10⋮ym0)F=(f1(t,y1,...,tm)⋮f1(t,y1,...,tm))

метод Эйлера:

{y′=F(t,y)y(t0)=y0

yi+1=yi+hF(ti,yi)

{amy(m)(t)+am−1y(m−1)+...+a1y′(t)+a0y(t)=f(t)y(t0)=y0⋮y(m−1)(t0)=y0(m−1)

y1(t)=y(t);y2(t)=y′(t)=;...;ym(t)=y(m−1)(t) - замены

{y'1=y2y'2=y3⋮y'm−1=ymy'm=f(t)−am−1ym−...−a0y1amy=(y1⋮ym)F=(y2⋮ymf(t)−am−1ym−...−a0y1am)

Устойчивость численных методов

[править]

1h∑j=0kαjyn+1−j=Φ(tn,yn+1−k,...,yn,yn+1,h)(1)

α0,α1,...,αk−1α0≠0

Чтобы облегчить задачу будем исследовать устойчивость по начальным данным (без правой части) y0,y1,...,yk−1

Определение. Метод устойчив по начальным данным, если

max0≤n≤N|yn∗−yn|≤k(τ)max0≤i≤k−1|εi|, где y∗ - решение (1) при возмущенных входных данных y0∗...yk∗. εi=yi∗−yi, i=0,k−1‾


Рассмотрим модельную задачу y′=0 на нуль-устойчивость

1h∑j=0kαjyn+1−j=0

(2)α0yn+1+α1yn+...+ykαn+1−k=0 - конечно-разностное уравнение. Стартовые точки yn+1−k...yn+1

Затем мы находим y0,y1,...,yk−1,yk,yk+1,...,yN

Возмущенная задача: y0∗...yk−1∗

α0yn+1∗+α1yn∗+...+αkyn+1−k∗=0

yi∗i=0,N‾ - проверим на устойчивость

α0εn+1+α1εn+...+αkεn+1−k=0

εn=ε‾qn, q некоторое число. ε‾ - некоторая первоначальная погрешность. Получим характеристическое уравнения k - шагового метода:

α0qk+α1qk−1+...+αk=0

P(q)=∑j=0kαkqk−j - характеристический полном Уравнение P(q)=0 имеет ровно k корней (кратных и простых). Пусть q1 - кратный корень кратности r, остальные простые:

q=c1(1)q1n+c1(2)nqn+...+c1(r−1)nr−1q1n+c2q2n+...+cmqmn(3)

Общее решение уравнения: |qi|<1,n→∞,|q|→0

Определение. Будем говорить, что выполнено корневое условие, если все корни полнинома устойчивости по модулю меньше 1, или на гранирце круга |qi|=1 кратных корней нет


Теорема. Если выполнено корневое условие, то метод является устойчивым (называется O-устойчивым)
Доказательство. из (3): |qi|<1εn=ε‾anεn→0 при n→∞

|qi|=1, но корень не кратен. εn=ε‾ - погрешность не растет.


Утверждение. Все методы Рунге-Кутты и методы Адамаса являются нуль устойчивыми методами.


Р-К: P(q)=q−1

А: α0qk+α1qk−1+0qk−2+...+0=0α0=1,α1=−1P(q)=qk−qk−1

Наличие O-устойчивости означает, что метод устойчив на конечном отрезке [t0,T]. Если T→∞ следует применять абсолютно устойчивые методы.

Модельная задача:

{y′=λyλ∈Cy(t0)=y0Reλ<0}⇒ Ассимптотическая устойчивость

1h∑j=0kαjyn+1−j=λ∑j=0kβj(λh)yn+1−jz=λh

∑j=0k(αj−zβj(z))yn+1−j=0γ=αj−zβj(z)

Для невозмущенной задачи: ∑j=0kγjyn+1−j=0y0,...,yk−1 - стартовые точки

P(q,z)=γ0qk+γ1qk−1+...+γk=0. Ищем q1,q2,...,qn (зависят от z)

Для явного метода Эйлера: yn+1−ynh=λyn

εn+1−(1+λh)εn=0q−(1+z)=0⇒q1=1+z|1+z|<1

−1<1+z<1⇒z∈(−2,0)z=λh;λ<0⇒0<h<2−λ

Метод называется устойчивым при ограничении на шаг h<2|λ|

Для неявного:

yn+1−ynh=λyn+1⇒11−λh=q⇒|1−λh|<1 - верно т.к. λ<0

Ограничений на h нет. Область устойчивости: |1−z|<1

Определение. Метод называется абсолютно устойчивым при конкретном значении z=λh если все его корни полинома устойчивости лежат внутри единичного круга, а на границе нет кратных корней. Множество точек z комплексной плоскости в которых выполнено корневое условие называется областью абсолютной устойчивости метода.


Определение. Численный метод называется A-устойчивым, если область абсолютой устойчивости включает в себя полу плоскость Rez<0


Неявный метод Эйлера A-устойчив

Понятие о жестких задачах

[править]
S=max⁡|Reλi(A)|min⁡|Reλi(A)|,S - число жесткости системы

s≫10 - система жесткая. Для таких задач надо применять A-жесткие методы.