Численные методы решения задачи Коши

на равномерной сетке (x₀ = a, x₁ , x₂ , … , x_m = b) отрезка [a, b] с шагом h = (b – a)/m являются методами Рунге-Кутта, если, начиная с данных (x₀ , y₀ ), решение ведется по следующим рекуррентным формулам:

(6)

Метод называют методом Рунге-Кутта порядка р, если он имеет р-й порядок точности по шагу h на сетке. Порядок точности р достигается с помощью формул (6) при определенных значениях коэффициентов c_j и d_j , j = 1, 2, … , p; коэффициент с₁ всегда полагают равным нулю. Эти коэффициенты вычисляются по следующей схеме:

1) точное решение (х₀ + h) и его приближение y₁ = y₀ + y₀(h) представляют в виде разложения по формуле Тейлора с центром х₀ вплоть до слагаемого порядка h ^p+1;

2) из равенств подобных членов при одинаковых степенях h в двух разложениях получают уравнения, решая которые находят коэффициенты c_j и d_j .

Метод Эйлера можно назвать методом Рунге-Кутта первого порядка. Действительно, для р = 1, с₁ = 0, d₁ = 1 формулы (6) преобразуются в соотношение (4):

,

или

Метод Рунге-Кутта второго порядка называют методом Эйлера-Коши, если р = 1, с₁ = 0, с₂ = 1, d₁ = d₂ = 1/2. Алгоритм метода Эйлера-Коши получается из формул (6):

. (7)

Для практической оценки погрешности решения можно применять правило Рунге, полагая в формуле (5) р = 2:

.

Метод Рунге-Кутта четвертого порядка называют классическим методом Рунге-Кутта, если р = 1, с₁ = 0, с₂ = с₃ = 1/2 , с₄ = 1, d₁ = d₄ = 1/6, d₂ = d₃ = 1/3. Из рекуррентных формул (6) получим алгоритм решения задачи Коши классическим методом Рунге-Кутта:

,

,

, (8)

,

.

Графиком приближенного решения является ломаная, последовательно соединяющая точки P_i(x_i , y_i ), i = 1, 2, …, m. С увеличением порядка численного метода звенья ломанной приближаются к ломаной, образованной хордами интегральной кривой y = (x), последовательно соединяющими точки (x_i , (x_i)) на интегральной кривой.

Правило Рунге (5) для практической оценки погрешности приближенного решения для численного метода четвертого порядка имеет вид:

.

Порядок выполнения лабораторной работы.

Пример. Найти решение задачи Коши для дифференциального уравнения первого порядка с начальным условиемна равномерной сетке отрезка[0; 1] один раз с шагом h = 0,2 , другой – с шагом 0,1 методами Эйлера, Эйлера-Коши и классическим методом Рунге-Кутта. Оценить погрешность численного решения по принципу Рунге. Сравнить численное решение с точным (x) = 2e ^x – x - 1.

Вид рабочего листа табличного процессора MS Excel приведен на рисунке.

1. В диапазоне A3:L10 оформляем вычисления для метода Эйлера. В диапазоне A3:А10 размещаем заголовки строк. Задаем узловые точки сетки: В3="0", С3="0.1", протягиваем диапазон В3:С3 до ячейки L3.

Программируем рабочий лист для вычислений с шагом h=0,2. Размещаем начальное условие : ячейкаВ5="1". Вычисляем угловой коэффициент в начальной точке f(x₀ , y₀ )=х₀+у₀: ячейка D4="=B3+B5". Вычисляем следующее приближенное значение функции у₁ = у₀ + h f(x₀ , y₀ ): ячейка D5="=B5+(D3-B3)*D4". Протягиваем диапазон D4:E5 в диапазон D4:L5 для вычислений в остальных узловых точках.

Программируем рабочий лист для вычислений с шагом h=0,1. Размещаем начальное условие : ячейкаВ7="1". Вычисляем угловой коэффициент в начальной точке f(x₀ , y₀ )=х₀+у₀: ячейка С6="=B3+B7". Вычисляем следующее приближенное значение функции у₁ = у₀ + h f(x₀ , y₀ ): ячейка С7="=B7+(C3-B3)*C6". Протягиваем диапазон С6:С7 в диапазон С6:L7 для вычислений в остальных узловых точках.

Вычисляем значение функции (x) = 2e ^x – x - 1 в узловых точках: ячейка В8="=2*EXP(B3)-B3-1" и протягиваем формулу в диапазон B8:L8. По формуле вычисляем оценку точности в узловых точках более крупной сетки: ячейка D9="=ABS(D7-D5)" и протягиваем диапазон D9:E9 до ячейки L9. По формуле вычисляем фактическое отклонение приближенного значения на сетке h/2 в узлах более крупной сетки: ячейка D10="=ABS(D8-D7)" и протягиваем диапазон D10:E10 до ячейки L10. Вычисления по методу Эйлера закончены.

2. В диапазоне A13:L30 оформляем вычисления для метода Эйлера-Коши. В диапазоне A13:A30 размещаем заголовки строк. Задаем узловые точки сетки: В13="0", С13="0.1", протягиваем диапазон В13:С13 до ячейки L13.

Программируем рабочий лист для вычислений с шагом h=0,2. Размещаем начальное условие : ячейкаВ20="1". Проводим расчет в 1-ой узловой точке. Вычисляем угловой коэффициент в начальной точке f(x₀ , y₀ )=х₀+у₀: ячейка D14="=B13+B20". Вычисляем коэффициент : ячейкаD15="=(D13-B13)*D14". Вычисляем аргументы коэффициента k₂: ячейка D16="=D13" (), ячейкаD17="=B20+D15" (y₀ + k₁ ). Вычисляем функцию f = х + у от полученных аргументов: ячейка D18 ="=D16+D17". Вычисляем коэффициент : ячейкаD19="=(D13-B13)*D18". Вычисляем следующее приближенное значение функции у₁ = у₀ + (k₁ + k₂ )/2 : ячейка D20="=B20+(D15+D19)/2". Протягиваем диапазон D14:E20 в диапазон D14:L20 для вычислений в остальных узловых точках.

Программируем рабочий лист для вычислений с шагом h=0,1. Размещаем начальное условие : ячейкаB27="1". Проводим расчет в 1-ой узловой точке. Вычисляем угловой коэффициент в начальной точке f(x₀ , y₀ )=х₀+у₀: ячейка C21="=B13+B27". Вычисляем коэффициент : ячейкаC22="=(C13-B13)*C21". Вычисляем аргументы коэффициента k₂: ячейка C23="=C13" (), ячейкаC24="=B27+C22" (y₀ + k₁ ). Вычисляем функцию f = х + у от полученных аргументов: ячейка C25 ="=C23+C24". Вычисляем коэффициент : ячейкаC26="=(C13-B13)*C25". Вычисляем следующее приближенное значение функции у₁ = у₀ + (k₁ + k₂ )/2 : ячейка C27="=B27+(C22+C26)/2". Протягиваем диапазон C21:C27 в диапазон C21:L27 для вычислений в остальных узловых точках.

Вычисляем значение функции (x) = 2e ^x – x - 1 в узловых точках: ячейка В28="=2*EXP(B13)-B13-1" и протягиваем формулу в диапазон B28:L28. По формуле вычисляем оценку точности в узловых точках более крупной сетки: ячейка D29="=ABS(D27-D20)/3" и протягиваем диапазон D29:E29 до ячейки L29. По формуле вычисляем фактическое отклонение приближенного значения на сетке h/2 в узлах более крупной сетки: ячейка D30="=ABS(D28-D27)" и протягиваем диапазон D30:E30 до ячейки L30. Вычисления по методу Эйлера-Коши закончены.

3. В диапазоне A33:L66 оформляем вычисления для классического метода Рунге-Кутта. В диапазоне A33:A66 размещаем заголовки строк. Задаем узловые точки сетки: В33="0", С33="0.1", протягиваем диапазон В33:С33 до ячейки L33.

Программируем рабочий лист для вычислений с шагом h=0,2. Размещаем начальное условие : ячейкаВ48="1". Проводим расчет в 1-ой узловой точке. Вычисляем угловой коэффициент в начальной точке f(x₀ , y₀ )=х₀+у₀: ячейка D34="=B33+B48". Вычисляем коэффициент : ячейкаD35="=(D33-B33)*D34". Вычисляем аргументы коэффициента k₂: ячейка D36="=B33+(D33-B33)/2" (), ячейкаD37="= B48+D35/2" (y₀ + k₁/2 ). Вычисляем функцию f = х + у от полученных аргументов: ячейка D38 = "=D36+D37". Вычисляем коэффициент : ячейкаD39="=(D33-B33)*D38". Вычисляем аргументы коэффициента k₃: ячейка D40="=B33+(D33-B33)/2" (), ячейкаD41="=B48+D39/2" (y₀ + k₂/2 ). Вычисляем функцию f = х + у от полученных аргументов: ячейка D42 = "=D40+D41". Вычисляем коэффициент : ячейкаD43="=(D33-B33)*D42". Вычисляем аргументы коэффициента k₄: ячейка D44="=D33" (), ячейкаD45="= B48+D43" (y₀ + k₃ ). Вычисляем функцию f = х + у от полученных аргументов: ячейка D46 = "=D44+D45". Вычисляем коэффициент : ячейкаD47="=(D33-B33)*D46". Вычисляем следующее приближенное значение функции у₁ = у₀ + (k₁ + 2k₂ + 2k₃ + k₄ )/6 : ячейка D48="=B48+(D35+2*D39+2*D43+D47)/6". Протягиваем диапазон D34:E48 в диапазон D34:L48 для вычислений в остальных узловых точках.

Программируем рабочий лист для вычислений с шагом h=0,1. Размещаем начальное условие : ячейкаВ63="1". Проводим расчет в 1-ой узловой точке. Вычисляем угловой коэффициент в начальной точке f(x₀ , y₀ )=х₀+у₀: ячейка C49="=B33+B63". Вычисляем коэффициент : ячейкаC50="=(C33-B33)*C49". Вычисляем аргументы коэффициента k₂: ячейка C51="=B33+(С33-B33)/2" (), ячейкаC52="=B63+C50/2" (y₀ + k₁/2 ). Вычисляем функцию f = х + у от полученных аргументов: ячейка C53="=C51+C52". Вычисляем коэффициент : ячейкаC54="=(C33-B33)*C53". Вычисляем аргументы коэффициента k₃: ячейка C55="=B33+(C33-B33)/2" (), ячейкаC56="=B63+C54/2" (y₀ + k₂/2 ). Вычисляем функцию f = х + у от полученных аргументов: ячейка C57 = "=C55+C56". Вычисляем коэффициент : ячейкаC58="=(C33-B33)*C57". Вычисляем аргументы коэффициента k₄: ячейка C59="=C33" (), ячейкаC60="=B63+C58" (y₀ + k₃ ). Вычисляем функцию f = х + у от полученных аргументов: ячейка C61="=C59+C60". Вычисляем коэффициент : ячейкаC62="=(C33-B33)*C61". Вычисляем следующее приближенное значение функции у₁ = у₀ + (k₁ + 2k₂ + 2k₃ + k₄ )/6 : ячейка C63="= B63+(C50+2*C54+2*C58+C62)/6". Протягиваем диапазон C49:C63 в диапазон C49:L63 для вычислений в остальных узловых точках.

Вычисляем значение функции (x) = 2e ^x – x - 1 в узловых точках: ячейка В64="=2*EXP(B33)-B33-1" и протягиваем формулу в диапазон B64:L64. По формуле вычисляем оценку точности в узловых точках более крупной сетки: ячейка D65="=ABS(D63-D48)/15" и протягиваем диапазон D65:E65 до ячейки L65. По формуле вычисляем фактическое отклонение приближенного значения на сетке h/2 в узлах более крупной сетки: ячейка D66="=ABS(D64-D63)и протягиваем диапазон D66:E66 до ячейки L66. Вычисления по классическому методу Рунге-Кутта закончены.