Как решить СДУ методом Рунге-Кутта?

Question

[Сергей Хлопов]

Здравствуйте, подскажите пожалуйста, делаю лабораторную работу, нужно было составить системы дифференциальных уравнений Колмогорова для определения вероятностей состояния системы:

И теперь нужно составить алгоритм программы для определения вероятностей состояний системы на отрезке времени [0;T] с шагом h, написано что можно решить это методом Рунге-Кутты.
И я так понимаю что бы решить данную систему нужно задать начальное условие для P:

Подскажите пожалуйста как решить данную систему ? В интернете смотрел но не могу разобраться с формулами которые там приведены, что-куда подставлять.

Заранее благодарю за ответ.

Comments

[lz961]

Даже в Википедии численная схема выписана в явном виде. Чего же вам ещё надо? Понять какой брать величину шага?

[Сергей Хлопов]

lz961,
Благодарю за комментарий.
А можете пожалуйста подсказать, сейчас пытаюсь разобраться, я так понимаю что для каждого уравнения которое у меня есть высчитываю: k1, k2, k3, k4.
x и y которые передаются в функцию в моём случае p и t если я верно понимаю.

Но не пойму как быть с "t", у меня в уравнениях есть к примеру p2(t) т.е. с изменение "t" должно и "p2" меняться ?
Но когда задавал начальное состояние системы p2 было 0.

По вот этим формулами делаю, это в интернете нашел, есть отличия от формул в Википедии

[lz961]

у вас игрики -- это p, икс это t. t у вас пробегает значения на сетке t0, t1=t0+h, t2=t1+h=t0+2*h, ...
если нужна вся история эволюции вероятностей, то храним её в массивах p1[N], p2[N], p3[N], p4[N], N = [T/h].
поскольку y теперь вектор из 4-х элементов, а время явно в уравнения не входит, делаем так:
k11 = h*(-.04*p1[i]+.04*p2[i])
k12 = h*(.01*p1[i] + .02*p3[i]+....)
k13 = ...
k14 = ...
k21 = h*(-.04*(p1[i]+k11/2) + .04*(p2[i]+k12/2))
....
p1[i+1] = p1[i] + 1/6*(k11+k21*2+k31*2+k41)
p2[i+1] = p2[i] + 1/6*(k12+k22*2+k32*2+k42)
....

[Сергей Хлопов]

lz961, Спасибо и подскажите ещё пожалуйста а почему в k21 идёт деление k11 на 2 ? Вроде h делим на 2

[lz961]

Сергей Хлопов, h относится к иксу, вместо которого у нас время, и он в явном виде при вычислении правой части уравнения не фигурирует.

[Сергей Хлопов]

lz961, Спасибо, и если можете подскажите ещё пожалуйста, немного запутался.
Когда составлял уравнение была такая вот формула:

Т.е. сумма всех полученных p должна получится 1, а теперь думаю, в любом же случае будет 1, даже наверно немного больше, если в одном из массивов p1, p2, p3, p4 у меня будет 1, т.к. начальное состояние системы задано, верно ?

[lz961]

Сергей Хлопов, в теории да, если правильно подобраны коэффициенты. На практике различие будет, и отличие суммы решений от 1 можно использовать для проверки корректности выбора шага h.

[Сергей Хлопов]

lz961, Понял, спасибо. И последний вопрос если не затруднит подскажите пожалуйста, я уже и так много вопросов задал)).
Вопрос по вот этому выражению:
k21 = h*(-.04*(p1[i]+k11/2) + .04*(p2[i]+k12/2))
во второй части (p2[i]+k12/2).

А в формуле же k1 делим на 2, т.е. по логике думал что вот так должно быть:
k21 = h * (0.04 * (p2[t] + k11 / 2) - 0.04 * (p1[t] + k11 / 2)); т.е. везде k1 для первого уравнения