Метод конечных разностей для решения алгоритма ОДУ - PullRequest
0 голосов
/ 06 апреля 2019

Я пытаюсь разработать алгоритм для метода конечных разностей, но я немного запутался.Рассматриваемый ODE: y '' - 5y '+ 10y = 10x, где y (0) = 0 и y (1) = 100.Поэтому мне нужен способ как-то получить коэффициенты, которые будут умножать "y_i" из отношения:

enter image description here

, а затем сохранить результирующие коэффициенты в матрице, которая будет матрицей системы, которую я решу через Гаусса-Иордана.Вопрос сводится к тому, как получить эти коэффициенты и переместить их в матрицу.Я думал о том, чтобы определить коэффициенты вручную, а затем просто ввести матрицу, но мне нужно сделать это для шагов размером 0,1, 0,001 и 0,001, так что это действительно нереальный вариант.

1 Ответ

2 голосов
/ 07 апреля 2019

Давайте предположим более общий случай ОДУ

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
Добро пожаловать на сайт PullRequest, где вы можете задавать вопросы и получать ответы от других членов сообщества.
...