график на фазовой плоскости
дано:
y(x0) = a
y'(x0) = b
найти:
y' = f(x,y)
решение:
x1 = x0 + d
y1 = y0 + d * y'0 = a + d * b
y'1 = (y1 - y0) / d
y2 = y1 + d * y'1 etc
недостаток - очень быстро копится систематическая ошибка, поэтому надежны только несколько первых значений полученной таким образом производной
дано: y(x0) = a0
y'(x0) = b0
найти: y' = f(x,y)
предположим что функция в окрестностях точки x0 разложена в ряд Тейлора:
f(x) = a0 + a1 * (x - x0) + a2 * (x - x0)^2 + a3 * (x - x0)^3 + ...
тогда
y'(x) = a1 + 2 * a2 * (x - x0) + 3 * a3 * (x - x0)^2 + ...
подставим в уравнение y'(x) = f(x) и приравняем соответствующие коэффициенты
например:
y' = y * x + 1
y(0) = 1
решение нельзя записать аналитически в виде зависимости y = f(x)
но представив y(x) в виде ряда Тейлора и подставив в дифуравнение:
a1 + 2 * a2 * x + 3 * a3 * x^2 + .. = 1 + a0 * x + a1 * x^2 + a2 * x^3 + ..
мы тут же получим что
a1 = 1
2 * a2 = a0 = y(0) = 1 откуда a2 = 1/2
3 * a3 = a2 откуда a3 = 1/6
...
можно получить сколь угодно точное решение в виде ряда. проверимся:
y = 1 + x + x^2/2 + x^3/6 + ...
y' = 1 + x + x^2/2 + ... = 1 + x * (1 + x + x^2/2 + x^3/6 + ..)
= 1 + x + x^2 + x^3/2 + O(x^4) все сходится при x << 1
t = 0, 1, 2, ...
s = ( x₁
x₂
.
xₙ )
s₁ = A * s₀
s₂ = A * s₁
..
sₙ - sk = A(n-k) * sk
старинное уравнение (операция) Ньютона
A * x = y ,
yi = xi+1 - xi
пусть есть последовательность временных квантов: t₀, ∆t + t₀, ...., k * ∆t + t₀
пусть есть вектор x = ( x₀
x₁
.
xₙ )
пусть yt = A * xt
тогда, полагая ∆t = 1, имеем
yt+1 - yt = ∆yt и матрица такого преобразования (полагая квант ∆t равным 1):
( y'₁ ( -1 1 ( x₀
y'₂ -1 1 x₁
y'₃ = -1 1 * A * x₂
. . .
y'ₙ ) -1 1 ) xₙ )
для разности второго порядка:
2 * yt - yt+1 - yt-1 = ∆∆yt
и матрица такого преобразования:
( y˝₁ ( 2 -1 ( x₀
y˝₂ -1 2 -1 x₁
y˝₃ = -1 2 -1 * A * x₂
.. .. ..
y˝ₙ ) -1 2 ) xₙ )
sn+1 - sₙ = A * s ₙ
матрица преобразования:
sn+1 = (I + A) * sₙ
в коде октав это будет что-то вроде того:
% forward Euler method
% u˝ + u = 0; u(0) = 1, u'(0) = 0
% u' = v
% v' = - u
% un+1 = un + ∆ * vn
% vn+1 = vn - ∆ * un
% |un+1| | 1 ∆| |un|
% | | = | | * | |
% |vn+1| |-∆ 1| |vn|
function fem (n)
x = y = zeros (n, 1) ;
d = 2 * pi / 100 ; % step
x(1) = 1 ; y(1) = 0 ; % initial conditions
A = [1 d ; -d 1] ;
for t = 2 : n
Z = A * [x(t-1) ; y(t-1)] ; x(t) = Z(1) ; y(t) = Z(2) ;
end
plot (x, y, 'linewidth', 2)
end
метод быстрый, но нестабильный, что легко видеть на графике - спираль раскручивается. сумма собственных значений равна трейсу матрицы (т.е. 2), а произведение равно детерминанту (т.е. больше единицы т.к. 1+d²>1). матрица несимметричная, значит собственные значения - комплексные и одно из них по модулю больше 1, откуда и следует нестабильность системы. таково се ля ви
разница между соседними витками спирали характеризует точность
ошибка на одном шаге алгоритма составляет u˝ * (∆t)² / 2
суммарная ошибка за N шагов будет
(1 + A * (∆t)²)N * err₁ + .. + (1 + A * (∆t)²)N-k * errk + .. + errₙ откуда следует критерий стабильности решения:
|1 + A * (∆t)²| ≤ 1
в общем случае (не для данного уравнения) если решение стабильно, то суммарная ошибка будет ограничена сверху значением:
u˝max * N * (∆t)²/2
заметив, что N * ∆t = T , где T - общее время наблюдения за системой (в тиках), получаем:
err ~ T/2 * ∆t * u˝max
sn+1 - sₙ = A * sn+1
ЗАМЕЧАНИЕ: метод неявный, потому что искомое есть как в lhs, так и в rhs, в отличие от прямого метода, где к моменту вычисления sn+1 значение sₙ было уже известно (либо из начальных условий на первом шаге, либо как результат предыдущего шага)
матрица преобразования:
sn+1 = (I - A)⁻ * sₙ
в октаве можно наваять что-либо в духе:
% backward Euler' method
% u˝ + u = 0 u(0) = 1 u'(0) = 0
% u' = v
% v' = - u
% un+1 - un = ∆ * vn+1
% vn+1 - vn = - ∆ * un+1
% |un+1 - un| | 0 ∆| |un+1|
% | | = | | * | |
% |vn+1 - vn| |-∆ 0| |vn+1|
function bem (n)
x = y = zeros (n, 1) ;
x(1) = 1 ; y(1) = 0 ; % initial conditions
d = 2 * pi / 100 ; % step
A = [0 d ; -d 0] ;
B = inv (eye (2, 2) - A) ;
for t = 1 : n
Z = B * [x(t) ; y(t)] ; x(t+1) = Z(1) ; y(t+1) = Z(2) ;
end
plot (x, y, 'linewidth', 2)
end
BEM стабильнее (по крайней мере в данном случае он сходится - спираль скручивается), чем FEM, но медленнее. сумма собственных значений матрицы равна нулю (трейс), а произведение собственных значений (детерминант) равен d², значит оба собственных значения - комплексные числа (сопряженные), а значит вещественная часть равна нулю, а модуль меньше 1 - поэтому колебания и затухают
стабильность метода определяется как |1 - 2 * a / Δt| ≤ 1
разница между соседними витками спирали характеризует точность
в отличие от метода Эйлера, который берет в рассмотрение Un или Un+1 данный метод рассматривает промежуточное значение на расстоянии Δt/2 как от un так и от un+1:
u' = a * u
( un+1 - un ) / Δt = a*un+1/2 + a*un/2
разобьем один шаг, на котором дифференциал не изменяет своего значения на два подшага
sn+1 - sₙ = A * sₙ
полагая, что следующее приращение значения такое же как и предыдущее:
sn+2 - sn+1 = sn+1 - Sₙ
получаем второе уравнение
sn+1 - sₙ = A * sn+1
теперь складываем левые и правые части и делим на два:
sn+1 - sₙ = 0.5 * A * (sn+1 + sₙ)
sn+1 - 0.5 * A * sn+1 = sₙ + 0.5 * A * sₙ
NB: метод неявный, потому что искомое есть как в lhs, так и в rhs
матрица преобразования:
sn+1 = (I - 0.5 * A)⁻ * (I + 0.5 * A) * sₙ
в октаве это может выглядеть как-то так:
% trapezoid method
% u˝ - u = 0
% u' = v
% v' = - u
% |u'| | 0 1| |u |
% | | = | | * | |
% |u˝| |-1 0| |u'|
%
% un+1 = (I - (∆t/2) * A)⁻ * (I + (∆t/2) * A) * un
function trapezoid (n)
x = y = zeros (n, 1) ;
x(1) = 1 ; y(1) = 0 ; % initial conditions
d = 2 * pi / 100 ; % step
A = [0 1 ; -1 0] ;
I = eye (2, 2) ;
B = (inv (I - (d / 2) * A)) * (I + (d / 2) * A) ;
for t = 2 : n
Z = B * [x(t-1) ; y(t-1)] ; x(t) = Z(1) ; y(t) = Z(2) ;
end
plot (x, y, 'linewidth', 2) ;
end
матрицу нужно вычислять лишь однажды, поэтому метод быстрый. транспонированная матрица равна обратной и значит все собственные числа лежат на единичной окружности комплексной плоскости. кроме того, трейс равен нулю и значит собственные числа чисто мнимые. именно по этой причине колебания не затухают и не взлетают
это же имеет и свои ограничения - все начинает работать не так в случае даже небольшой нелинейности задачи
стабильность метода определяется как |[1 + a*Δt/2]/[1 - a*Δt/2]| ≤ 1. для отрицательных значений а метод абсолютно стабильный, для положительных - нестабильный
точность метода на каждом шаге: err ~ |u'''| * (∆t)³ / 12
при решении системы дифуров одна часть решается используя значения из предыдущего шага, а вторая часть решается уже используя это вновь полученное значение
пусть u˝ + u = 0
положим, что u' = v
тогда v' = -u
получившаяся система:
un+1 = uₙ + ∆t * vₙ
vn+1 = vₙ - ∆t * un+1
или
un+1 = uₙ + ∆t * vₙ
∆t * un+1 + vn+1 = vₙ
( 1 0 ( un+1 ( 1 ∆t ( uₙ
∆t 1 ) * vn+1 ) = 0 1 ) * vₙ )
преобразуем. найдем обратную для левой матрицы:
1 0 | 1 0
0 1 | - ∆t 1
умножим получившуюся на правую и левую части:
( un+1 ( 1 ∆t ( uₙ
vn+1 ) = -∆t 1 - ∆t*∆t ) * vₙ )
можно немножко покодить:
% leap-frog method
% u'' + u = 0 u(0) = 1 u'(0) = 0
% u' = v
% v' = - u
% un+1 = un + ∆vn
% vn+1 = vn - ∆un+1 = vn - ∆(un + ∆vn)
% |un+1| | 1 ∆ | |un|
% | | = | | * | |
% |vn+1| | -∆ 1-∆*∆ | |vn|
function leap (n)
x = y = zeros (n, 1) ;
x (1) = 1 ; y (1) = 0 ; % initial conditions
d = 2 * pi / 100 ; % step
A = [1 d ; -d 1-d*d] ;
for t = 2 : n
Z = A * [x(t-1) ; y(t-1)] ; x(t) = Z(1) ; y(t) = Z(2) ;
end
plot (x, y, 'linewidth', 2) ;
end
un+1 - uₙ un+1 - 2 * uₙ + un-1
------------ + ------------------------ = A * un+1
Δt 2 * Δt
на первом шаге надо получить значение в следующей точке используя FEM. дальше - по алгоритму
практически в технических расчетах хватает достаточно точного и стабильного метода Рунге-Кутта второго порядка
суть метода заключается в том, что раскладывают значение функции в ряд Тейлора: значение функции в нулевой точке, первая и вторая (взвешенная) производные в сумме дают достаточно точный обобщенный результат
набросок кода выглядит приблизительно так:
% runge-kutta second-order method
% err ~ C * ∆^3
% ∆ = A * x + A * A * x * 0.5
% u˝ = - u
% u' = v
% v' = - u
% |u'| | 0 1| |un|
% | | = | | * | |
% |v'| |-1 0| |vn|
function runge_kutta (n)
x = y = zeros (n , 1) ;
x (1) = 1 ; y (1) = 0 ; % initial conditions
d = 2 * pi / 100 ; % step
A = [0 d ; -d 0] ;
for t = 2 : n
Z0 = [x(t-1) ; y(t-1)] ; % the function value in x0
Z1 = A * Z0 ; % the first derivation in x0
Z2 = A * Z1 ; % the second derivation in x0
Z = Z0 + Z1 + 0.5 * Z2 ; % f(x) = f(x0) + f`(x0) * d + f˝(x0) * d²/2
x(t) = Z(1) ; y(t) = Z(2) ;
end
plot (x, y, 'linewidth', 2) ;
end
результат применения этого кода:
* * *
на небольших матрицах все считается достаточно быстро. если матрицы большие, то их сначала надо преобразовать (диагонализировать, к примеру) и потом считать уже с оптимальными матрицами. но это уже другая история
посчитать диффузию (без конвекции) двумя способами - конечных разностей и среднего значения:
уравнение диффузии: Δ u = f * du/dt
listing octave
то же самое - в хаскеле
listing haskell