Давайте предположим более общий случай ОДУ
c1 * y''(x) + c2 * y'(x) + c3 * y(x) + c4 * x = 0
с граничными условиями
y(0) = lo
y(1) = hi
И вы хотите решить это для x ∈ [0, 1]
с размером шага:h = 1 / n
(где n + 1
- количество образцов).Мы хотим решить для yi = y(h * i)
.Диапазон yi
от i ∈ [0, n]
.Для этого мы хотим решить линейную систему
A y = b
Каждая внутренняя часть yi
будет налагать одно линейное ограничение.Следовательно, у нас есть n - 1
строк в столбцах A
и n - 1
, соответствующих неизвестным yi
.
Чтобы настроить A
и b
, мы можем просто скользить окном над нашимнеизвестно yi
(я предполагаю, что индексация начинается с нуля).
A = 0 //the zero matrix
b = 0 //the zero vector
for i from 1 to n - 1
//we are going to create the constraint for yi and store it in row i-1
//coefficient for yi+1
coeff = c1 / h^2 + c2 / h
if i + 1 < n
A(i - 1, i) = coeff
else
b(i - 1) -= coeff * hi //we already know yi+1
//coefficient for yi
coeff = -2 * c1 / h^2 - c2 / h + c3
A(i - 1, i - 1) = coeff
//coefficient for yi-1
coeff = c1 / h^2
if i - 1 > 0
A(i - 1, i - 2) = coeff
else
b(i - 1) -= coeff * lo //we already know yi-1
//coefficient for x
b(i - 1) -= c4 * i * h
next