Численное моделирование эволюционных задач.

Преподаватель: Бобков Владимер Васильевич

Какая - то задача, кот.связана с процессом эволюции, отслеживанием процесса эволюции

Главный Вопрос : как математической эволюционной задаче поставить в соответствие численную модификацию - вычислетельные методы
	предметный области не рассматриваются.

эволюционная задача
	нет строгово определения
	1. Задача, которая описывает процесс, развивающийся во времени
	2. ДУ, которое содержит только первые производные по времени.
	очень мало условий - только одно начальное условие и есть много условий
	типичные примеры:
	задача Коши для системы ОДУ

u'[x]=f[x,u[x]], u[t]=y; x>t (t=0).

Чем неприятно? Если мы хотим отследить эволюцию x, то мало дополнительных условий - в отличие от краевой задачи "второй конец" свободен, не закреплён.
Эта "простая" задача для численного моделирования более сложная. То, что важно в ДУ нужно не потерять на этапе приближения.
Сделали шаг - то, что до этого забыли, продолжаем. Но если чуть-чуть отклонились - стали на другое решение. Если нет асимпотической устойчисовсти, то ошибка не исправляется, а накапливается.
Нужно так провести замену исходной задачи, чтобы наиболее важное не пропало. А что это важное? - большой вопрос.
	Если бы были граничные условия - была бы проверка в конце. Но у нас только начальное.

Пример: u'[x]=f[x]  - нет зависимости от х
ошибки накапливаются. Если решаем на ограниченном отрезке, то модуль итоговой ошибки будет не больше суммы модулей отдельных ошибок.
Отдельные ошибки: (((τ)^2)/2)u''[ξ]
Сделав N шагов: ((T-(Subscript[t, 0]))/N)=τ;  Если устремить τ к нулю, то и итоговая ошибка устремиться к нулю.
Если считать не на конечном промежутке, то ошибка будет накапливаться всё равно - даже если разностная схема второго порядка.

Пример:

u'[x] = A u[x] + a; u[t]
(u^(стационарное))=-(A^(-1))a

Если мы не выйдем на бесконечности на -(A^(-1)), мы принципиально исказим вычисление. Есть методы сколь угодного порядка точности, которые не выводят на бесконечности на правильный результат. Они годятся если исследовать ограниченные отрезки, то, может, и нормально. Но если исследовать процесс на бесконечности, то нужны хорошые методы.

Точное решение:
u[t+τ]=(e^(Aτ))u[t]+(∫_0^τ)(e^(Aξ))dξ a

Рассмотрим приближённый вариант: два члена от экспоненты (e^Az)=I+((A z)/(1!))+(((A^2)(z^2))/(2!)) - явный метод, метод ломаных.
	(OverHat[y])=(I+A τ)y + (∫_0^τ)(I+A ξ)dξ a = (I+A τ)y + (τ+((τ^2)/2)A)a = E y +S a

(y^(стационарное))=(y^(стационарное))+τ A (y^(стационарное))+τ a + ((τ^2)/2)A τ; (y^(стационарное)) = -(A^(-1))a - (τ/2)a

((-τ)/2)a - ошибка. Даже при малых шагах τ её можно сделать большой при помощи вектора свободных членов a.

Схемы, которые, казалось бы, нормальны могут привести к плачевному результату.

Эту схему исправить легко, но реальные могут быть сложны.

Традиционно извесные вещи

Метод переменных направлений, например, не соблюдает порядок "выгорания"  гармоник (сначала большее x в (e^(-x)), потом мешьние и т.д.)

Даже неявные схемы сохраняют предельное положение только тогда, когда рассчинаты _точно _ (а их точно не считают --  по 2 - 3 итерации).

Фундамент

Есть исходная эволюционная задача. Предполагаем, что решение сущестует и единственно.
Задача: создать численную модель, пригоную для анализа этой задачи.
Превый вопрос: вопрос о том, какие источники ошибок и какие составляющие ошибок приближённого решения.

Основные источники и составные части погрешностей и способы их математического описания.

Есть точное решение, и есть приближённое решение, полученное с "отягчающими обстоятельствами".
Погрешность - разность между точным и приближённым решениями.
Пусть есть одно уравнение.
u'[x] = f[x, u[x]], u[t] = y; (* Система ОДУ. Решение существует и единственно. *)
Численный метод:
 (OverHat[y])=y  +τ p[t,τ,y,(OverHat[y])]; (* p - приближение к производной *)
Погрешности
1. Неустранимая - та, которая зависит только от погрешности начальных условий (несовершенство датчиков)
	Неустранимые ошибки уже не удалишь. Они следуют законам дифференциального уравнения.
	(* если df/du<0, то траектории сближаются; =0 - параллельны; >0 - расходятся *)
	Всё может быть описано на языке неустранимой погрешности.
	Неустранимыми также считают погрешности в других входных данных (помимо начального условия)
	Каждый шаг допускается некоторая ошибка. На следующем шаге она рассматривается как неустранимая.
2. Вычислительная - погрешность округления. Разрядная сетка ограничена.
	будем уделять им мало внимания
	нельзя делить на числа, близкие к нулю
	нельзя вычитать числа близкие , чтобы главное значение не исчезло и точность не уменьшилась.
	нельзя проводить рекурсивные процедуры, умножая на числа больше 1 (или соответствующие матрицы)
(* вероятностные методы широко используются в оценке вычислительных методов. Воеводин. *)
3. погрешность замены исходной задачи - наша основная тема. Ту задачу, которая нам не по силам заменяют берущейся  - основной источник ошибок.
 Та часть погрешности, которая зависит только от апроксимации (замены) называют погрешность метода.
 Чем же оценить эту погрешность (замены)?
Limit[p[t,τ,u,u[t+τ]]  , τ->0] = f[t, u] для любого решения u (* необходимое условие *)
Подставим точное решение в модель.: (OverHat[u])=u+τ p[t,τ,u,(OverHat[u])]+r; r - невязка.
	это опосредованное описание точности решения
	хороша она тем, что можно выяснить её порядок малости
		 r = ±((τ^2)/2)u''[ξ]
	Если r -> 0 (при τ -> 0) стремится со скоростью O[(τ^(q+1))], то решение будет стремиться с O[(τ^q)]
	Если метод явный, то r[t,r]=ϵ[t+τ]=u[t+τ]-(OverHat[y]) - невязка равна неточности решение.
		маленькая невязка - маленькая ошибка
	недостатки
		нужно знать ξ и u''[t]
		для неявного метода невязка не отражает точность решения
Конкретный пример:
Пусть есть набор уравнений:
1. u'[x]=100u[x]-10(x^2)+0.2x-1 (* неустойчивость *)
2. u'[x]=-100u[x]+10(x^2)+0.2x+1 (* асимптотическая устойчивость *)
3. u'[x]=0.2x (* устойчивое (но не асимпотически) *)
Начальные условия: u[0]=0.01 (* для всех *)
Точное решение (для всех): u[x]=0.1(x^2)+0.01
Метод: неявный метод Эйлера: (OverHat[y])=y+τ (OverHat[f])
Невязка: r=-((τ^2)/2)u''[ξ]=-0.1(τ^2)
Пусть τ=0.01. Невязка: -0.00001.
Локальная ошибка: u[0.01]-(OverHat[y])->+∞ (* даже на первом шаге неточность большая (хотя невязка маленькая) *)
	(OverHat[y])=y+τ 100(OverHat[y])-…; (OverHat[y])=(…/(1-100τ)); (* для первого уравнения *)
Для τ<0.01 будет более-менее OK. А для τ>0.01 будет как во втором, устойчивом уравнении => _существенное_ искажение (на самом деле расходится, а по модели сходится). А по невязке всё выглядит нормально.

Другой, альтернативный метод оценивания замены.

u'[x] = f[x, u[x]], u[t] = y;
 (OverHat[y])=y  +τ p[t,τ,y,(OverHat[y])];

Задача решается на отрезке t<=x<=t+τ
Введём понятие производной разностного метода.
	Свяжем с нашей моделью функцию непрерывного аргумента:
	y[x]=y+(x-t)p[t,(x-t),y,y[x]]=y+(x-t)q[x,y[x]]

	Продифференцируем её:
	 y'[x] = q[x, 
   y[x]] + (x - t) ((Subscript[q, x])[x, y[x]] + (Subscript[q, y])[x, y[x]] y'[x])=ϕ[x,y[x]]

Поставим в соответствие нашему исследуемому методу дифференциальное уравнение: 
	y'[x]=ϕ[x,y[x]]
	
	Это похоже на наше исходное уравнение u'[x]=f[x, u[x]]

y'[x]=(((1-(x-t)(Subscript[q, y])[x,y[x]]))^(-1))(q[x,y[x]]+(x-t)(Subscript[q, x])[x.y[x]])

Применим к той нашей задаче:
y'[x]=(((1-(x-t)100))^(-1))(…) будет расходимость (когда x=0.01)
наше решение непрерывно, а тут уже при t=0.01 разрыв производной.
	
Можно говорить о невязке....

y'[x] = ϕ[x, y[x]] - приближённое дифференциальное уравнение.
u'[x] = f[x, u[x]] - точное дифференциальное уранение

r[x]=f[x,u[x]]-ϕ[x,u[x]] - невязка точного ....
ρ[x]=f[x,y[x]]-ϕ[x,y[x]] - невязка приближённого решения. Она может быть просто вычислена.

Опасности остаются. Но теперь они видны (нули в знаменателях)
ρ[x] - главная часть производной локальной погрешности. (∫_t^(t+τ))ρ[c]dc

Пример:Возьмём задачу
u'[x]=A u[x]. Пусть A - симметричная, чтобы был базис и собственные вектора.
(OverHat[u])=Exp[E τ]u; u=(∑_(i=1)^n)(Subscript[C, i])(ξ^i)
(OverHat[u])=(∑_(i=1)^n)(Subscript[C, i])(ξ^i)Exp[(Subscript[λ, i])τ]
(OverHat[y])=y+τ A (OverHat[y])
(OverHat[y])=(((I-τ A))^(-1))y
y=u
(OverHat[y])=(∑_(i=1)^n)(Subscript[C, i])(ξ^i)(1/(1-τ (Subscript[λ, i])))
спектральная функция,множитель перехода E=(1/(1-λ τ))
(рисунок) E=E(λ,τ)  Для λ<-(2/τ) уже плохо.
Теперь для неявного:
E=((1)/(1+τ λ))

Plot[{{(*явный метод*)1+λ τ},(*неявной метод*){(1/(1-λ τ))},{Exp[λ τ]}}/.τ->0.5,{λ,-4,4}](*график переходного множителя в зависимости от λ*)

Для явного метода при λ< -(2/τ) прямая ниже экспоненты на графике проблемы
для неявного метода там хорошо - убывает почти как экспонента
но для больших лямбда при некотором τоно неверно, потом ешё и ложная устойчивость может быть
y'[x]=ϕ[x,y[x]]; y[t]=y
Как сделать чтобы была исключена большая погрешность при маленькой невязке?
Наша погрешность: ϵ'[x]=u'[x]-y'[x]=f[x,u[x]]-ϕ[x,y[x]]=
=(({{f[x,u[x]]-ϕ[x,u[x]]}, {==r[x]}}))+(({{ϕ[x,u[x]]}, {==y[x]+ϵ[x]}}))-ϕ[x,y[x]]
e[t]==0; 
ϵ'[x]=ρ[x]+f[x,((u[x])/(y[x]+e[x]))]-f[x,y[x]];

u'[x]=f[x,u[x]],u[t]=y;
Если явный метод
	(OverHat[y])=y+τ p [t,t,y];
(OverHat[u])=y+τ p [t,t,y]+r;
(OverHat[u])-(OverHat[y])=r[t,τ]; (*для неявного *)
Для неявного: ......=ϕ[x,(*без y[x]*)];r[x]=ϵ'[x];ϵ[x]=(∫_t^x)r[z]dz
r мы не можем вычислить - оно на точном решении. Но ρ[x]≈ϵ'[x].
ϵ[x]≈(∫_t^x)ρ[z]dz - главная часть ошибки
Добавка к ϵ'[x]: f[x,((u[x])/(y[x]+e[x]))]-f[x,y[x]];
Pros: 1. Больше информации, 2. Чувствительно к "ляпам" с пвесдоустойчивостью.

пока что мы только более грамотно опывасали пределы методов....

<offtopic>Пример: уравнение A u==f. Есть единственное решение.
Если A разложить на ((Subscript[A, 1])+(Subscript[A, 2 ])+...+(Subscript[A, r])) .
 ((Subscript[A, 1])+(Subscript[A, 2 ])+...+(Subscript[A, r])) u==f. Пока всё нормально
Но когда берут  (Subscript[A, 1])(Subscript[u, 1])+(Subscript[A, 2 ])(Subscript[u, 2])+...+(Subscript[A, r]) (Subscript[u, r])==f - уже fail. (оно делалось чисто в расчёте на малый шаг. И слишком малым он может показаться даже для машины. Делать акцент на малость шага преподаватель не рекомендует. С бесконечно малым шагом бесконечно далеко не уйдёшь.</offtopic>

Представление в виде последовательности задач

Нельзя и пойти задом наперёд? Описывая ... некоторой дифференциальной задачей. Сделать апроксимирующие задачи, решения которых приближаются к нашей, и которые решать легко (численно).
Как мы проведём на шаге апроксимацию?
Нужно предложить такие задачи вида y'[x]=ϕ[x,y[x]]
1. легко и надёжно решать
2. можно построить последовательность задач чтобы решения сходились к решению нашей задачи (на шаг).
Проблема не в том, что много считать, а что нужно надёжно считать.
невязка: ρ[x]=f[x,y[x]]-y'[x]; 
неточность: ϵ'[x]=ρ[x]+f[x,y[x]+ϵ[x]]-f[x,y[x]]; ϵ[t]=0.
Простейший пример:
(Subscript[u', i+1])[x]=f[x,(Subscript[u, i])[x]]=(Subscript[a, i+1])[x]; (Subscript[u, i+1])[t]=y; (Subscript[u, 0])=y. Если процесс сойдётся, то будет то, что нужно. Если процес сходится, то z'[x]=f[x,z[x]].
(Subscript[u, i+1])[x]=y+(∫_t^x)f[z,u'[z]]dz; Получилась теорема Пикара. Условие сходимости извесно.
Почему он не используется? 1. Медленно сходится, 2. Вычислительная сложность
u'[x]=A u[x]; A=Aᵀ<0; (OverHat[u])=Exp[A τ] u;
(Subscript[u, p])[x]=(∑_(i=0)^p)((((x-t))^i)/(i!))(A^i)y=(Subscript[E, p])[A,x-t]y;
Спектральные значения: (Subscript[E, p])[(Subscript[λ, i]),τ]=(∑_(i=0)^p)((τ^i)/(i!))(Subscript[λ, i])
y=(∑_(i=1)^n)(Subscript[d, i])(ξ^i); (OverHat[u])=(∑_(i=1)^n)(Subscript[d, i])(ξ^i)Exp[(Subscript[λ, i])τ]==(∑_(i=1)^n)(Subscript[d, i])(γ^i)[({{t+τ}, {== x}})]
(Subscript[u, p])[x]=(∑_(i=1)^n)(Subscript[d, i])(ξ^i)(Subscript[E, p])[(Subscript[λ, i]),x-t]=(∑_(i=1)^n)(Subscript[d, i])(Subscript[q, i])[x];  q - гармоники
Решение - линейная комбинация гармоник матрицы A.

Чтобы соблюжадлся порядок выгорания нармоник, 
0<Exp[(Subscript[λ, i])(x-t)]<Exp[(Subscript[λ, j])(x-t)]<1
0  < (Subscript[E, p])[(Subscript[λ, i]), x - t] < (Subscript[E, p])[(Subscript[λ, j]),x-t]<1

(Subscript[E, p ])- первые p+1 членок ряда Макларена для экспоненты.

Если в таком методе взять слишком большое τ, то (Subscript[E, p]) будет стремиться к +∞ или -∞ (в зависимости от чётности p) (при p=∞ всё сходится чётко).
А если шаг τ слишком маленький, то из-за недостаточной машинной точности младшие гармоники "зависнут", не будут выгорать.
(* все гармоники у нас в решении в сумме - неточность будет в обоих случаях *)
Из-за этого процесс Пикара не очень используется.
(* усиленная асимптотическая монотонность - спектральная монотонность (не только нельзя выгорать медленнее самой медленной, но и для всех остальных должен быть такая же очерёдность выгорания)*)

Простым штрихом улучшаем метод Пикара.

Цена: отказ от независимости (Subscript[u', i+1])[x]от u'[x].

u'[x]=f[x,u[x]]; ставиться в соответствие y'[x]=ϕ[x,y[x]].
Последовательность ϕ[x]
(Subscript[u, i+1])'=f[x,(Subscript[u, i])[x]]; (Subscript[u, 0])[x]=y; (Subscript[u, i+1 ])[x]=y+(∫_t^x)f[z,(Subscript[u, i])[z]]dz - это обычный процесс. Но плохо: нарушаются требования: нужно брать мелкий шаг - плохо. Ещё один недостаток: на кажом шаге всё решение строится заново. Хотелось бы получить прибл. решение, а потом поправки-поправки (всё уменьшающиеся).

Модифициуем чтобы можно было уточнять.
Будем считать, что y'[x]=ϕ[x,y[x]] - неизвестная присутствует и в правой части.
Для начала выберем другое начальное приближение: (Subscript[u, 0])'[x]=A (Subscript[u, 0])[x]+(Subscript[a, 0]). Как выбрать A и (Subscript[a, 0])? Разложим по Тейлору:f[x,u[x]]=f[t,y]+(Subscript[f, x])[t,y](x-t)+(Subscript[f, u])[t,y](u[x]-y)+(Subscript[δ, 2]), Norm[(Subscript[δ, 2])] ∈ O[(((x-t))^2)].
A=(Subscript[f, u])[t,y]; (Subscript[a, 0])=f-A y  (* (Subscript[f, x])[t,x] лучше не включать в (Subscript[a, 0]) - теряется порядок малоти *)

(Subscript[u, 0])[x] - решение этого уровнения:
(Subscript[u, 0])[x]=Exp[A(x-t)]y+(∫_t^x)Exp[A(x-z)](Subscript[a, 0])dz=Exp[A(x-t)]y+(∫_0^(x-t))Exp[A ξ]dξ (Subscript[a, 0])=
=Exp[A(x-t)]y + (Exp[A(x-t)] - I)(A^(-1))(Subscript[a, 0])

Какова погрешность такой замены? (Subscript[ϵ, 0])'[x]=f[x,u[x]]-A (Subscript[u, 0])[x]-(Subscript[a, 0])=(Subscript[ρ, 0])[x]+(Subscript[g, 0])[x];
(Subscript[ρ, 0])[x]=f[x,(Subscript[u, 0])[x]]-A (Subscript[u, 0])[x]-(Subscript[a, 0]); (Subscript[g, 0])[x]=f[x,u[x]]-f[x,(Subscript[u, 0])[x]];
(Subscript[ρ, 0])[t]=0; (Subscript[g, 0])[t]=0. Е ещё (Subscript[ϵ, 0])[t]=0. И ещё получили (Subscript[ϵ, 0])'[t]=0.
Значит (Subscript[ϵ, 0])[x]=(Subscript[ϵ, 0])[t]+(x-t)(Subscript[ϵ, 0])'[t]+(Subscript[δ, 2])[x-t] = (Subscript[δ, 2])[x-t]. - получился метод решения первого порядка точности.

Он допускает реализацию, метод явный, СЛАУ решать не нужно. Нужно чтобы было ((A)^(-1)) и чтобы удобно было матричную экспоненту вычислять.

Далее: мы не будет по-обычному приближаться к решению, а как-то по-другому уточнять (решая (Subscript[u, 0])'[x]=A (Subscript[u, 0])[x]+(Subscript[a, 0])).
Считаем, что (Subscript[u, 0])[x]уже имеем, будем искать поправки.

(Subscript[ϵ, o])'[x]=(Subscript[ρ, 0])[x]+(Subscript[g, 0])[x];  (Subscript[ϵ, 0])[x]=({{(∫_t^x)(Subscript[ρ, 0])[z]dz}, {== (Subscript[δ, 0])[x]}})+({{(∫_t^x)(Subscript[g, 0])[z]dz}, {== (Subscript[ϵ, 1])[x]}})

(Subscript[ϵ, 0])[x]== (({{(Subscript[δ, 0])[x]}, {==u[x]-(Subscript[u, 0])[x]}}))+(({{(Subscript[ϵ, 1])[x]}, {==u[x]-(Subscript[u, 1])[x]}})); (Subscript[u, 1])[x]=(Subscript[u, 0])[x]+(Subscript[δ, 0])[x]

Norm[(Subscript[ϵ, 1])[x]]∈O[(((x-t))^3)] - наше решение уже второго порядка точности

(Subscript[ϵ, 1])'[x]== u'[x]-(Subscript[u, 1])'[x]== f[x,u[x]]-(({{(Subscript[u, 1])'[x]}, {==f[x,(Subscript[u, 0])[x]]}}))==
==(({{(Subscript[ρ, 1])[x]}, {== f[x,(Subscript[u, 1])[x]]-(({{(Subscript[u, 1])'[x]}, {== (Subscript[u, 0])'[x]+(Subscript[δ, 0])'[x] }}))}}))+(({{(Subscript[g, 1])[x]}, {== f[x,u[x]]-f[x,(Subscript[u, 1])[x]]}}))

(Subscript[u, 1])'[x]== (Subscript[u, 0])'[x]+(Subscript[δ, 0])'[x]==(Subscript[u, 0])'[x]+f[x,(Subscript[u, 0])[x]]-(Subscript[u, 0])'[x]

Потом(Subscript[ϵ, 1])'[x]=(Subscript[ρ, 1])[x]+(Subscript[g, 1])[x];  (Subscript[ϵ, 1])[x]== (({{(∫_t^x)(Subscript[ρ, 1])[z]dz}, {== (Subscript[δ, 1])[x]}}))+(({{(∫_t^x)(Subscript[g, 1])[z]dz}, {== (Subscript[ϵ, 2])[x]}})) - третий порядок точности и т.д.

Опишем сам процесс приближений :

u'[x]=A (Subscript[u, 0])[x]+(Subscript[a, 0]) - сначала

Потом (Subscript[u, i+1])[x]=(Subscript[u, i])[x]+(Subscript[δ, i])[x]
(Subscript[δ, i])[x]=(∫_t^x)(Subscript[ρ, i])[z]dz
(Subscript[ρ, i])[x]=f[x,(Subscript[u, i])[x]]-f[x,(Subscript[u, i-1])[x]];
(Subscript[u, i+1])'[x]=f[x,(Subscript[u, i])[x]];

Преимущества :
1. Не нужно решать СНЛАУ - многие законы  теряются если приближённо
2. Нет тех проблем, которые мы описывали
3. Сходимость : если сходится, то туда, куда нам нужно. Правила когда сходится определены.

Основная проблемы : 
  1. хороший способ матричной экспоненты
2. хорошая решалка СТЛОД.

...


Предположения:
1. Гладкость исходного уравнения
2. Теория .... .... специализированные методы их решения

1. Эта задача - специфический неявный метод с первым порядком точности.
2. Сходимость как у метода Пикара

Замечание : если встретится  f[x,u[x]]=A u[x]+(({{a[x]}, {переменное,а не постоянное}})), то имея точное решение (Subscript[u, 0])[x] для постоянного случая (с (Subscript[a, 0])) можно без бесконечного приближения, в один шаг перейти к решению переменного :
(Subscript[ϵ, o])'[x]=A (Subscript[ϵ, 0])[x]+a[x]-(Subscript[a, 0]); (Subscript[ϵ, 0])[t]=0
(Subscript[ϵ, 0])[x]==(∫_t^x)Exp[A(x-z)](a[z]-(Subscript[a, 0]))dz== (Subscript[u, 0])[x]+∫...
(опять упирается в проблему матричной экспоненты)

Пусть мы сделали начальный шаг: нашли (Subscript[u, 0])[x]
Для (Subscript[ϵ, 0])[x] имеем (Subscript[ρ, 0])[x]+(Subscript[g, 0])[x].
Попробуем как находить (Subscript[ϵ, 0])[x] не из простого, а из дифференциального уравнения.
(Subscript[ϵ, 0])[x]== (Subscript[ρ, 0])[x]+(Subscript[g, 0])[x]==A (Subscript[ϵ, 0])[x]+(Subscript[ρ, 0])[x]+(({{(Subscript[g, 0])[x]-A (Subscript[ϵ, 0])[x]}, {== (Subscript[q, 0])[x]}}))

(Subscript[q, 0])[x]==(Subscript[g, 0])[x]-A (Subscript[ϵ, 0])[x] = f[x,(({{u[x]}, {==(Subscript[u, 0])+(Subscript[ϵ, 0])}}))]-f[x,(Subscript[u, 0])[x]]-A (Subscript[ϵ, 0])[x]=(({{(Subscript[f, u])}, {матрица Якоби}})[x,(Subscript[u, 0])[x]]-A)(Subscript[ϵ, 0])[x]+(({{(Subscript[s, 2])[x]}, {norm∈O[((Norm[(Subscript[ϵ, 0])[x]])^2)]}}))

Norm[(Subscript[ϵ, 0])[x]] ∈O[(((x-t))^2)]
Norm[(Subscript[g, 0])[x]]∈O[(((x-t))^2)]
Norm[(Subscript[q, 0])[x]]∈O[(((x-t))^3)]

(Subscript[ϵ, 0])[x]+(({{(∫_t^x)Exp[A(x-z)](Subscript[ρ, 0])[z]dz}, {==(Subscript[δ, 0])[x]}}))+(({{(∫_t^x)Exp[A(x-z)](Subscript[q, 0])[z]dz}, {==(Subscript[ϵ, 1])[x]; norm∈O[(…^4)]}}))

Итак, теперь на первом приближении у нас уже 2 порядка точности

(Subscript[ϵ, 1])'[x]==A (Subscript[e, 1])[x]+(Subscript[q, 0])[x]. 
А хотелось бы A (Subscript[e, 1])[x]+(Subscript[ρ, 1])[x]+(Subscript[q, 1])[x]; где ({{(Subscript[ρ, 1])[x]=f[x,u[x]]-(Subscript[u, 1])'[x]}, {(Subscript[q, 1])[z]=f[x,u[x]]-f[x,u[x]]-A ϵ[x]}})

Нужно доказать: (Subscript[ρ, 1])[x]+(Subscript[q, 1])[x]==q[x]

(Subscript[ρ, 1])[x]+(Subscript[q, 1])[x]==f[x,u[x]]-(Subscript[u, 1])[x]+({{f[x,u[x]]}, {==u'[x]}})-f[x,(Subscript[u, 1])[x]]-A (Subscript[ϵ, 1])[x]=(Subscript[ϵ, 1])'[x]-A ϵ[x]=(Subscript[q, 0])[x]

Norm[(Subscript[q, 1])[x]]∈O[(((x-t))^(5 ))]; (Subscript[ϵ, 1])[x]==({{(∫_t^x)Exp[A(x-z)](Subscript[ρ, 1])[x]dz}, {==(Subscript[δ, 1])[x]}})+({{(∫_t^x)Exp[A(x-z)](Subscript[q, 1])[z]dz}, {==(Subscript[ϵ, 2])[x]}})

На каждом рекурсионном приближении получаем + 2 порядка точности.

Опишем этот уточнённый процесс дифференциальных привязок :

u'[x]=A (Subscript[u, 0])[x]+(Subscript[a, 0]) - сначала

Потом (Subscript[u, i+1])[x]=(Subscript[u, i])[x]+(Subscript[δ, i])[x]
(Subscript[δ, i])[x]=(∫_t^x)Exp[A (x-z)](Subscript[ρ, i])[z]dz (*а было просто (Subscript[ρ, i])[z] *)
(Subscript[ρ, i])[x]=f[x,(Subscript[u, i])[x]]-(Subscript[u, i])'[x];(Subscript[u, i])'[x]==(Subscript[u, i-1])'[x]+(Subscript[δ, i-1])'[x]==(Subscript[u, i-1 ])'[x]+(Subscript[ρ, i-1])[x]+A (Subscript[δ, i-1])[x]==(Subscript[u, i-1 ])'[x]+f[x,(Subscript[u, i-1])[x]]-(Subscript[u, i-1 ])'[x]+A (Subscript[δ, i-1])[x]=A (Subscript[u, i])[x]+f[x,(Subscript[u, i-1])[x]]-A (Subscript[u, i-1])[x]==A (Subscript[u, i])[x]+(Subscript[a, i])[x]; 

(Subscript[a, i])[x]==f[x,(Subscript[u, i-1])[x]]-A (Subscript[u, i-1])[x]==(({{f[x,(Subscript[u, i-1])[x]]-(Subscript[u, i-1])'[x]}, {==(Subscript[ρ, i-1])[x]}}))+(Subscript[u, i-1])'[x]-A (Subscript[u, i-1])[x]==
==(*опустимся на одну рекурсию ниже *) == (Subscript[ρ, i-1])[x]+A (Subscript[u, i-1])[x]+(Subscript[a, i-1])[x]-A (Subscript[u, i-1])[x]==(Subscript[ρ, i-1])[x]+(Subscript[a, i-1])[x]
(Subscript[ρ, i])[x]=(Subscript[a, i+1])[x]-(Subscript[a, i])[x]=f[x,(Subscript[u, i])[x]]-f[x,(Subscript[u, i-1])[x]]-A (Subscript[u, i-1])[x]-A (Subscript[u, i])[x]+A (Subscript[u, i-1])[x]==...
({{(Subscript[u, i])'[x]}, {z'[x]}})=A (Subscript[u, i])[x]+f[x,(Subscript[u, i-1])[x]]-A (Subscript[u, i-1])[x]=f[x,({{(Subscript[u, i-1])[x]}, {z[x]}})]+({{A (Subscript[δ, i-1])[xs]}, {0}}).

Основное преимущество : точность повышается на 2 порядка на каждом шаге (а не на один), при этом поправки δ ... .. и вместо ρ, которое всё равно нужно считать просто экспонента.

Проблемы и требования к явным и неявным методам

(* телефон преподавателя: 2027684 *)

Характерные проблемы традиционных численных методов.
На примере явного и неявного Эйлера

u'[x] = A u[x] + a; u[t] = y
Законы/требования:
1. единственное положение равновесия: (если невырождена) (u^(ст))=-(A^(-1))a.

u[x]==-(A^(-1))a + Exp[A x] c; (OverHat[u])== Exp[a τ]y+((∫_t^τ)Exp[A ξ]dξ a == (Exp[A τ]-I)(A^(-1))a)
(* дожно работать даже когда вырожденная матрица *)0<E[(Subscript[λ, i]),τ]<1 для (Subscript[λ, i])<0, E[(Subscript[λ, i]),τ]>1 для (Subscript[λ, i])>0

2. Спектральная монотонность: Если (Subscript[λ, i])<(Subscript[λ, j]), то 0<E[(Subscript[λ, i]),τ]<E[(Subscript[λ, j]),τ]

Show[(Plot[#[[1]]/.{τ->0.5},{λ,-5,5},PlotRange->{-4,6},PlotStyle->#[[2]]]&/@(({{Exp[λ τ], Red}, {1+λ τ , Green (*explicit Euler*)}, {(1/(1-λ τ)), Blue (*implicit Euler*)}, {-1, Black}, {1, Black}})))]

(* absciss: lambda. Goal: to be between 0 and 1 (for λ<0), prefferably near the red line. Wrong sign - epic fail.*) (* ордината: множитель перехода *)

3'. Проблема подвисания медленных  гармоник (уже рассматривалось: слишком малый шаг, малая лямбда => произведение (1+(Subscript[λ, i])τ) ниже машинной точности. Вместо мизерного изменения получаем вообле отсутствие изменений).

Проблема жёсткости : в неявных лучше. Если точность большая - проблема подвисания. Точность маленькая - расходимости.

4. Ещё одна проблема неявного эйлера: λ<0, вроде как-бы всё хорошо. На графике гипербола почти приближает экспоненту. Рассмотрим ошибку ϵ[(Subscript[λ, i]),τ]=Exp[(Subscript[λ, i])τ]-1/(1-(Subscript[λ, i])τ). Ошибка, вроде, маленькая.  Если (Subscript[λ, i])τ==0, то ошибка вроде равна нулю. Если τ->∞, то ϵ тоже ->0. Но по пути экстремум. При этом (Subscript[λτ, э])≈-2.5,(Subscript[ϵ, э])≈-0.2

Plot[(-Exp[(Subscript[λ, i])τ]+1/(1-(Subscript[λ, i])τ)/.(Subscript[λ, i])->#)&/@{-0.5,-3,-10},{τ,0.1,5}]

(*Абсцисса: τ. Уменьшая τ мы не можем быть уверены, что мы уменьшаем ошибку. Выбор τ вообще ограничен в неявном методе *)

Потому что E[(Subscript[λ, i]),τ]->0 при τ->∞,(Subscript[λ, i])<0 - требование L-устойчивости.
Но это есть недостаток L-устойчивых методов.

... .. Почему бы не отказаться от лишней устойчивости (в том числе опасной ложной устойчивости) для получения чего-то полезного.

Рассмотрим ещё проблему неоднородности. В конце же нужно иметь -(A^(-1))a.

(OverHat[y]) ==E[A,τ] y+g[A,τ,a]
(OverHat[y])== -E (A^(-1))a+(∑_(i=1)^n)(Subscript[C, i])(ξ^i)Exp[(Subscript[λ, i])t]E[(Subscript[λ, i]),τ]+g=(({{g-E (A^(-1))a}, {хотелось бы чтобы это было -(A^(-1))a}}))+(({{......}, {уже на стационаре, ==0}}))

Утверждение : если стою на стационаре, то после шага нужно чтобы на стационаре и остался.

g==(E-I)(A^(-1))a

(OverHat[y]) ==E[A,τ] y+g[A,τ,a]=y+Q y+g;

g==(E-I)(A^(-1))a=>-(Q^(-1))g==-(A^(-1))a.

(y^(ст))=(y^(ст))+Q (y^(ст))+g; (y^(ст))=-(Q^(-1))g;

(∫_0^τ)Exp[A ξ]dξ a == (Subscript[g, точное])
g=(SuperStar[S])a; (SuperStar[S])≈(∫_0^τ)Exp[A ξ]dξ; (SuperStar[S])=(∫_0^τ)(SuperStar[E])[A,ξ]dξ; (SuperStar[S])a = (E-I)(A^(-1))a; E = I+(SuperStar[S])A. Это условие достаточно для сохранения положения равновесия равновесия.

(SuperStar[S])=(∫_0^τ)(SuperStar[E])[A,ξ]dξ; (SuperStar[E])≈Exp[A ξ];

Кроме сохражения положения равновесия ещё будет выполняться свойство дифференциальной согласованности :
(u')[x]=A(u')[x]; y'[x]=(SuperStar[E])[A,x-t](A y+a)y'[t]; y[x]=E[A,x-t]y+(SuperStar[S])[A,x-t]a

Реализация этого метода

Излишняя устойчивость пропадёт, появится спектральная монотонность,
Хотим монотонность по E, нужно следать монотонность и по (SuperStar[E]), так как как выбрал (SuperStar[E]), так и ..........
Exp[A z]≈I+((z A)/(1!))+...+(((z^p) (A^p))/(p!))
(SuperStar[(Subscript[E, m])])[λ,ξ(*шаг*)]== (∑_(i=0)^m)({{(Subscript[a, i])}, {(*это число*)}})(((1+(λ/μ)))^i); μ≥Norm[A]. Правее μ не должно быть собственных значений.

(* там пишут в спектральном виде. Метод: (((1+...))^i)*)

(((δ^k)(SuperStar[(Subscript[E, m])])[-μ,ξ])/((δλ)^k))==(ξ^k)Exp[-μ ξ] для k ∈ {0,…,m}.

(Subscript[a, 0])==Exp[-μ ξ]; (Subscript[a, 1])=((μ ξ)/(1!))Exp[-μ ξ];...;(Subscript[a, i])=((((μ ξ))^i)/(i!))Exp[-μ ξ]

0<(SuperStar[(Subscript[E, i])])≤(SuperStar[(Subscript[E, i+1])])≤Exp[λ ξ] - требование монотонности

(((δ^k)(SuperStar[(Subscript[E, m])]))/(δ (λ^k)))=(ξ^k).........

Show[Plot[#[[1]],{λ,-4,4},PlotStyle->#[[2]]]&/@(({{Exp[λ 0.4], Red}, {Exp[-3*0.4] (* (SuperStar[(Subscript[E, 0])]) *), Black}, {Exp[-3*0.4](1+0.4(λ+3))(* (SuperStar[(Subscript[E, 1])]) *), Black}, {Exp[-3*0.4](1+0.4(λ+3)+(((0.4^2)(((λ+3))^2))/(2!)))(* (SuperStar[(Subscript[E, 2])]) *), Black}, {Exp[-3*0.4](1+0.4(λ+3)+(((0.4^2)(((λ+3))^2))/(2!))+(((0.4^3)(((λ+3))^3))/(3!)))(* (SuperStar[(Subscript[E, 3])]) *), Black}})),Graphics[Circle[{-3(* -μ *),Exp[-3*0.4]},0.1]]]

(* Ось абсцисс λ. то, что левее -μ нас не интересует. То, что дальше делаем всё ближеприближает экспоненту *)

(Subscript[E, m+1])=I+(A/μ)(∑_(i=0)^m)(Subscript[a, i,1])(((I+(A/μ)))^i)

(OverHat[y])=(Subscript[E, m+1])y+(SuperStar[(Subscript[S, m])])a

(*при любом τ, поднимая m, можно получить сколь угодно точное решение. И наоборот, завфиксировав m, можно уменьшать τ. Полезные свойства будут соблюдаться. *)

(*Вопросы, литература - потом будут*)