Блок-схема программы вычисления решения дифференциального уравнения по методу Рунге-Кутта.

Начало

Описание Описание

процедуры- процедуры-

-функции Ввод х₀, у₀, -функции

rk(x₀,y₀,h,m) h₀, E, b y′ =f(x, y)

x₀, y₀ – начальные условия

h0 – шаг таблицы

n = ] (b- x₀)/ h₀[ + 1 E погрешность

b – конец отрезка

интегрирования [a, b]

i = 1, n n – число точек в таблице

решений

] [ - целая часть числа

h = h₀

m = 1

y = rk() Вызов процедуры-функции

y = rk(x₀,y₀, h,m)

y1 = y, h = h/2

x = x₀, y = y₀, m = 2m

Вызов процедуры-функции

y = rk() y = rk(x₀,y₀, h,m)

да |y-y1|>15E/16

нет

Вывод

x₀, y₀, h, m

x₀ = x₀ + h₀, y₀ = y

end

Блок-схема процедуры-функции метода Рунге-Кутта.

rk(x, y, h, m) x, y – начальная точка

h – шаг интегрирования

m – число дроблений шага

j = 1, m

Вычисление у

по формулам

Рунге-Кутта

x = x + h

rk = y Значение решения ОДУ

в очередной точке.

Расчетные формулы для метода Рунге- Кутта.

П р и м е р .Дано дифференциальное уравнение y′ = x – y c начальным условием y(0) = 1.5. Методом Рунге-Кутта с точностью до Е =0. 1 составить таблицу значений частного решения этого уравнения на интервале [0, 1.5]. Примем шаг таблицы, равным h₀ = 0.25.

Программа этой задачи на языке QBASIC имеет вид

DEF FNur (x, y) = x - y

DEF FNrk (x, y, h, m)

FOR j = 1 TO m

k1 = FNur(x, y)

k2 = FNur(x + h / 2, y + h * k1 / 2)

k3 = FNur(x + h / 2, y + h * k2 / 2)

k4 = FNur(x + h, y + h * k3)

y = y + (k1 + 2 * k2 + 2 * k3 + k4) * h / 6

x = x + h

NEXT j

FNrk = y

END DEF

x0 = 0

y0 = 1.5

h0 = .25

E = .01

b = 1.5

n = INT((b - x0) / h0) + 1

PRINT x0, y0, h0, n

FOR i = 1 TO n

h = h0

m = 1

y = FNrk(x0, y0, h0, m)

m1: y1 = y

h = h / 2

x = x0

y = y0

m = 2 * m

y = FNrk(x, y, h, m)

IF ABS(y - y1) > 15 * E / 16 THEN GOTO m1

WRITE "xi=", x0, "yi=", y0, "h=", h, "m=", m, "i=", i

x0 = x0 + h0

y0 = y

NEXT i

END

Результаты работы программы

x y h m (число дроблений шага)

0 1.5 0.125 2

0.25 1.1970 0.125 2

0.5 1.0163 0.125 2

0.75 0.9309 0.125 2

1.00 0.9196 0.125 2

1.25 0.9663 0.125 2

1.50 1.0578 0.125 2

В таблице взяты два резервных десятичных знака.

Решение методом Рунге-Кутта можно осуществить, используя пакет MathCAD.

Начальное значение y = y₀ , это обозначается ORIGIN =1. Далее обращаетесь к встроенной функции f(x) в разделе «Решение дифференциальных уравнений», функция «rkfixed»

В первой позиции вставляем у, во второй и третьей - границы изменения х, в четвертой – число точек деления интервала, в пятой позиции - функцию D(x,y), являющуюся правой частью дифференциального уравнения. Далее – «Z=», появляется таблица решений дифференциального уравнения. Первый столбец – значения переменной х, второй – значения у

Можно построить график решения. Вторая линия – график производной.