Донецкий национальный технический университет
Опубликован: 15.03.2007 | Доступ: свободный | Студентов: 5968 / 2125 | Оценка: 4.11 / 3.78 | Длительность: 12:32:00
Специальности: Математик
Лекция 12:

Компьютерное моделирование и решение нелинейных уравнений

< Лекция 11 || Лекция 12: 12345678910

Метод Рунге - Кутта 2-го порядка (модифицированный метод Эйлера)

Отбросим в (12.4) члены ряда, содержащие h3, h4, h5:.

Тогда

y(x_i+h)=y(x_i)+h\cdot y'(x_i)+\frac{h^2}{2!} \cdot y"(x_i). ( 12.6)

Чтобы сохранить член ряда, содержащий h2, надо определить вторую производную y"(xi).Ее можно аппроксимировать разделенной разностью 2-го порядка

y"=(x_i)=\frac{\Delta y'}{\Delta x}= \frac{y'(x_i+h)-y'(x_i)}{h}

Подставляя это выражение в (12.6), получим

y=(x_i+h) =y(x_i) + h\cdot y'(x_i) + \frac{h}{2} \cdot \frac{y'(x_i+h)-y'(x_i)}{h}=\\
=y(x_i) + \frac{h}{2}y'(x_i) + \frac{h}{2}y'(x_i+h).

Окончательно, модифицированная или уточненная формула Эйлера имеет вид:

y_{i+1}=y_i + \frac{h}{2}f(x_i,y_i) + \frac{h}{2}f(x_{i+1},y_{i+1}). ( 12.7)

Как видно, для определения функции y(x) в точке i+1 необходимо знать значение правой части дифференциального уравнения f(xi+1, yi+1) в этой точке, для определения которой необходимо знать предварительное значение yi+1.

Для определения предварительного значения yi+1 воспользуемся формулой Эйлера. Тогда все вычисления на каждом шаге по модифицированной или уточненной формуле Эйлера будем выполнять в два этапа:

На первом этапе вычисляем предварительное значение y_{i+1}^Э по формуле Эйлера

y_{i+1}^Э = y_i + hf(x_i,y_i).

На втором этапе уточняем значение y=i+1 по модифицированной или уточненной формуле Эйлера

y_{i+1} = y_i + \frac{h}{2}f(x_i,y_i) + \frac{h}{2}f(x_{i+1},y_{i+1}^Э).

Точность метода определяется отброшенными членами ряда Тейлора (12.4), т.е. точность уточненного или модифицированного метода Эйлера на каждом шаге \approx h^3.

Рассмотрим геометрический смысл модифицированного метода Эйлера.

Так как

f(x_i,y_i) = y'(x) = \tg\alpha_i,\\
f(x_{i+1},y_{i+1}) = y'(x_{i+1}) = \tg\alpha_{i+1},
то модифицированную формулу Эйлера можно представить в виде:

y_{i+1} = y_I + \frac{h}{2}\tg\alpha_i + \frac{h}{2}\tg\alpha_{i+1},

где

\tg \alpha_i - тангенс угла наклона касательной к искомой функции у(х) в начальной точке каждого шага,

\tg \alpha_{i+1} - тангенс угла наклона касательной к искомой функции у(х) в конечной точке каждого шага.

Геометрический смысл модифицированного  метода Эйлера

Рис. 12.12. Геометрический смысл модифицированного метода Эйлера

Здесь

P1 - накопленная ошибка в (i+1)й точке по методу Эйлера,

P2 - накопленная ошибка в (i+1)й точке по модифицированному методу Эйлера.

Как видно из рис.12.11, в первой половине каждого шага, то есть на участке [xi, xi+h/2], искомая функция y(x) аппроксимируется прямой, которая выходит из точки (xi, yi) под углом, тангенс которого \tg \alpha_i = f(x_i,y_i).

Во второй половине этого же шага, т.е. на участке [xi + h/2,xi + h], искомая функция y(x) аппроксимируется прямой, которая выходит из точки с координатами

x=x_i + h/2,\\
y=y_i = \frac{h}{2} \cdot f(x_i,y_i)

под углом, тангенс которого

\tg \alpha_{i+1} = f(x_{i+1},y_{i+1}^Э).

В результате в модифицированном методе Эйлера функция у(х) на каждом шаге аппроксимируется не одной прямой, а двумя.

Алгоритм модифицированного метода Эйлера можно построить в виде двух программных модулей: основной программы и подпрограммы МELER, реализующей метод (рис. 12.13).

Схема алгоритма модифицированного метода Эйлера

Рис. 12.13. Схема алгоритма модифицированного метода Эйлера

Здесь

(x,y) -при вводе начальная точка, далее текущие значения табличной функции,

h -шаг интегрирования дифференциального уравнения,

b -конец интервала интегрирования.

< Лекция 11 || Лекция 12: 12345678910
Равиль Султанов
Равиль Султанов

В уравнениях движения кривошипно-шатунного механизма вместо обозначения радиуса кривошипа "r" ошибочно записан символ "γ" (гамма).

P.S. Может быть это слишком очевидно, но не упомянуто, что угол поворота кривошипа φ считается малым.

Александр Никитин
Александр Никитин

Добрый день.

В расчете параметра Т4 xi суммируется с величиной h/2 ?