Добавил:
Опубликованный материал нарушает ваши авторские права? Сообщите нам.
Вуз: Предмет: Файл:

Пособие

.pdf
Скачиваний:
127
Добавлен:
28.01.2022
Размер:
2.37 Mб
Скачать

Погрешность метода Эйлера связана с величиной шага интегрирования отношением e1 =C1h2, где C1– произвольная постоянная.

Пример 6.2-1. Решить методом Эйлера ОДУ y = 2x/y с начальными

условиями x0 = 1и y0 = 1на отрезке [1;1.4]с шагом h = 0.2.

x

 

1.2;

y

 

y

 

h y

1

0.2

2 1

1.4;

 

 

0

 

 

1

 

1

 

0

 

 

1

 

 

 

 

 

 

 

 

 

 

 

 

 

 

x

 

1.4;

y

 

y

h y 1.4 0.2

1 1.2

1.742.

2

2

 

 

 

 

 

 

1

1

 

 

 

1.4

 

 

 

 

 

 

 

 

 

 

 

 

 

6.3. Методы Рунге-Кутты

Методы Рунге-Кутты – это группа методов, широко применяемых на практике для решения ОДУ. В этих методах при вычислении значения искомой функции в очередной точке хi+1 используется информация о предыдущей точке хi, yi. Методы различаются объемом вычислений и точностью результата.

Порядок метода Рунге-Кутты определяется кратностью вычисления значения производной искомой функции f(x,y) на каждом шаге. В соответствии с этим метод Эйлера является методом Рунге-Кутты 1-го порядка, поскольку для получения очередного значения yi+1 функция f(x,y) вычисляется один раз в предыдущей точке хi, yi. В методах Рунге-Кутты более высоких порядков для вычисления очередного значения искомой функции в точке хi+1 значение правой части уравнения y’= f(x,y) вычисляется несколько раз. Таким образом, порядок метода Рунге-Кутты определяется количеством вычисляемых производных, требуемых для получения очередной точки решения ОДУ.

Метод Рунге-Кутты 2-го порядка (усовершенствованный метод Эйлера). Вычисление значения искомой функции в точке хi+1 проводится в два

этапа. Сначала вычисляют вспомогательную величину yi 1 по методу Эйлера:

y

i 1

y

h f(x

,y

).

 

i

i

i

 

(6.3-1)

Затем значение производной искомой функции в точке (xi+1,yi+1) используется для вычисления окончательного значения функции:

 

 

 

 

 

 

 

 

 

 

y

 

y

 

h

f(x i ,yi ) f(xi 1,yi 1)

.

(6.3-2)

i 1

i

 

 

 

2

 

 

 

 

 

 

 

 

 

 

 

 

Подставляя (6.3-1) в (6.3-2), окончательно получим расчетную формулу метода Рунге-Кутты 2-го порядка:

51

 

 

y

 

h

 

 

 

) f(x

 

h,y

 

hf(x ,y

 

y

i 1

i

 

f(x

,y

i

i

)) ,

 

 

2

 

i

i

 

 

i

i

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

где

i 0,1,2,...,n 1.

 

 

 

 

(6.3-3)

Этот метод также называют методом прогноза и коррекций.

находят грубое приближение

yi 1

по методу Эйлера (прогноз),

уточненное значение yi+1(коррекция).

В общем виде формулу (6.3-3) можно представить как

Сначала а затем

y

 

y

 

 

h

(k

 

k

 

 

),

 

i 1

i

2

1

2

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

где

 

 

 

k

1

f(x

,y

),

 

 

 

 

 

 

 

 

 

 

i

 

 

 

i

 

 

 

 

k

2

f(x

i

h,y

i

hk

).

 

 

 

 

 

 

 

 

 

 

 

 

1

 

(6.3-4)

Метод Рунге-Кутты второго порядка имеет наглядную геометрическую интерпретацию (рис.6.3-1). Построение проводится следующим образом: yi 1

определяется пересечением перпендикуляра, восстановленного из точки xi+1 c касательной L1, проведенной к кривой y(x) в предыдущей точке i,yi). Затем в

точке

(x

i 1

, y

i 1

)

проводится прямая

L2 с тангенсом угла наклона, равным

 

 

 

f (xi 1, yi 1 )

 

 

 

(xi 1, yi 1 ) под углом, тангенс

. Прямую L проводят через точку

которого находим усреднением значений тангенсов углов наклона прямых L1 и

L2. Прямая L проводится параллельно

L через точку i,yi). Ее пересечение с

перпендикуляром, восстановленным из

точки хi+1, и дает уточненное значение

yi+1.

 

 

 

 

 

 

 

 

 

 

Рис. 6.3-1

Погрешность метода Рунге-Кутты второго порядка связана с величиной шага интегрирования отношением e2 =C2h3, где C2– произвольная постоянная.

52

Пример 6.3-1. Решить методом Рунге-Кутты второго порядка ОДУ y = 2x/y с начальными условиями x0 = 1 и y0 = 1 на отрезке [1;1.4] и шагом h = 0.2.

1.x1 x0 h 1 0.2 1.2,

y

 

y

 

h f(x

 

,y

 

) 1 0.2

2 1

1.4,

 

 

1

0

 

0

 

 

 

 

 

 

 

0

 

 

 

 

 

1

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

y

y

 

 

h

(f(x

 

 

) f(x ,y

))

1 0.1(2 1

2 1.2

) 1.3714.

0

 

0

 

 

1

 

 

2

 

 

 

1

1

 

 

 

1.4

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

2.x2 x1 h 1.2 0.2 1.4,

y

 

y

 

h f(x ,y

 

) 1.3714 0.2

2 1.2

1.7214,

 

 

 

 

2

1

 

 

 

 

 

 

 

 

 

 

 

1

1

 

 

 

 

 

 

1.3714

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

y

 

y

 

 

h

(f(x ,y

 

 

) f(x

,y

 

)) 1.3714 0.1(

2 1.2

 

2 1.4

) 1.709061

1.7091.

2

1

 

 

 

2

 

 

 

 

 

 

2

1

1

 

2

 

 

 

 

1.3714

 

1.7214

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

x

0

1,

 

 

 

y

0

1,

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

x

1

1.2,

 

 

y

1

1.3714,

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

x

2

1.4,

 

 

y

2

1.7091.

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Обобщая вышесказанное, можно записать общую формулу для решения ОДУ первого порядка методами Рунге-Кутты:

y

i 1

y

h (x

,y

,h),

 

i

i

i

 

(6.3-5)

где Ф – линейная функция аргументов x, y, h и f(x,y), которая может быть представлена как

(xi,yi,h) P1k1 P2k2 ... Pnkn,

где

k1

f(xi,yi );

 

 

k2

f(xi α2h,yi

α2hk1);

(6.3-6)

k3

f(xi α3h,yi α3hk2 );

 

.........................................

 

kn

f(xi αnh,yi

αnhkn 1).

 

Величина n в (6.3-6) определяется порядком метода, а коэффициентам2, 3, … , n, Р1, Р2, … ,Pn подбирают такие значения, которые обеспечивают минимальную погрешность. Так, для метода Рунге-Кутты четвертого порядка (n=4) получена расчетная формула при следующих коэффициентах:

2= 3=1/2, 4=1, P1 = P4=1/6, P2 = P3 =2/6.

Подставив значения коэффициентов в (6.3-4), имеем

53

y

i 1

 

где

y

 

 

h

(k

 

 

 

i

6

1

 

 

 

 

 

 

 

 

 

 

 

 

k

1

f(x

,

 

 

 

 

 

i

 

 

k

2

 

f(x

i

 

 

 

 

 

 

 

 

k

3

 

f(x

i

 

 

 

 

 

 

 

k

4

 

f(x

i

 

 

 

 

 

 

 

2k

2

2k

3

k

4

),

 

 

 

 

 

y

);

 

 

 

 

 

 

i

 

 

 

 

 

 

 

h / 2,yi hk1 / 2);

h / 2,yi hk2 / 2);

h,yi hk3 ).

(6.3-7)

Геометрическая интерпретация этого метода очень сложна и потому не приводится.

Погрешность метода Рунге-Кутты четвертого порядка значительно

меньше методов первого и второго порядков и пропорциональна величине h (e4 =C4h5).

Пример 6.3-2. Решить методом Рунге-Кутты четвертого порядка ОДУ y = 2x/y с начальными условиями x0 = 1 и y0 = 1 на отрезке [1;1.4] с

шагом h = 0,2.

x1 1.2; y1 1 0.26 (2 2 1.833 2 1.859 1.7495) 1.3711;

x2 1.4; y2 1.3711 0.26 (1.7518 2 1.6862 2 1.69 1.69393) 1.7089.

Сведем в таблицу результаты решения уравнения y =2x/y методами Рунге-Кутты, соответственно, первого (y1i), второго (y2i) и четвертого (y4i) порядков и сравним с результатами, полученными точным методом (yi).

Таблица 6.3-1

 

хi

 

 

y1i

 

 

y2i

 

 

y4i

 

 

yi

 

 

 

 

 

 

 

 

 

 

 

 

1

 

 

1

 

 

1

 

 

1

 

 

1

 

 

1.2

 

 

1.4

 

 

1.3714

 

 

1.37115

 

 

1.37113

 

 

1.4

 

 

1.74286

 

 

1.7091

 

 

1.7089

 

 

1.7088

 

На практике для обеспечения требуемой точности (при использовании любого приближенного метода решения ОДУ) применяется автоматический выбор шага методом двойного просчета. При этом в каждой точке хi по

формуле, соответствующей выбранному методу, производится расчет yiс шагом h (yi(h)) и с шагом h/2(yi(h/2)). Цель двойного просчета состоит в том, чтобы для

каждой точки численного решения эти значения отличались на величину, не превышающую заданной погрешности . В этом случае общая формула для оценки погрешности решения ОДУ методами Рунге-Кутты имеет следующий вид:

54

y

 

(h)

y

(h / 2)

 

 

 

 

 

 

 

i

 

 

 

i

ε,

 

 

2

 

1

 

 

 

p

 

 

 

(6.3-8)

где p – порядок метода Рунге-Кутты. Эта формула называется также

правилом Рунге.

Если | yi(h))- yi(h/2)|< , то шаг для следующей точки выбирается равным h, иначе шаг уменьшается вдвое и продолжается уточнение yi в точке хi.

6.4. Решение ОДУ n-го порядка

Методы, рассмотренные выше, позволяют найти численное решение ОДУ только первого порядка. Однако они применимы и к уравнениям n-го порядка. Для этого ОДУ n-го порядка предварительно приводится к системе n уравнений первого порядка.

Пусть, например,

требуется решить ОДУ второго порядка

y '' f(x,y,y ') , с

начальными условиями

 

0

,

 

0

 

0

,

y (x

0

0 .

 

 

x x

 

 

y(x

 

) y

 

 

 

) y

 

Обозначим z=y’. В результате подстановки в исходное уравнение получим систему двух уравнений первого порядка

 

y z;

 

 

z f(x,y,z)

с двумя неизвестными функциями

y(x

0

) y

0

,

z(x

0

) z

0

y

 

 

 

 

0 .

,

y y(x)

и z z(x)и начальными условиями

В общем виде система уравнений может быть представлена в виде

y f (x,y,z);

 

1

 

z f

(x,y,z).

 

2

 

(6.4-1)

Решением системы (6.4-1) являются две функции y(x) и z(x) , из которых y(x) - решение исходного уравнения второго порядка. Выбрав, например, метод Эйлера, приближенное решение системы (6.4-1) можно найти с помощью двух рекуррентных формул:

y

i 1

y

i

h f (x

,y

,z

);

 

 

 

 

1

i

i

i

 

z

 

z

 

 

h f

(x

,y

,z

).

 

i 1

i

 

 

2

i

i

i

 

Пример 6.4-1. Дано обыкновенное дифференциальное уравнение

второго порядка y'' 2xy y' при

начальных условиях x0 0 ,

y0 1,

 

2

y0

на отрезке [0;0.4] с шагом h 0.2 .

55

Обозначим

 

, тогда

z y

системы ОДУ первого порядка

ОДУ второго порядка можно записать в виде

z 2xy;

y z

с начальными условиями

x

0

 

0

,

y

0

1

 

,

z

0

 

2

.

Применим метод Эйлера для решения системы ОДУ

x

0.2;

1

 

x

2

0.4;

 

 

и т.д.

zi 1 zi hf1(xi,yi,zi )yi 1 yi hf2(xi,yi,zi )

z

 

 

z

 

h(2x

 

y

 

 

z

 

)2 0.2(2 0 1 2) 1.6;

 

1

 

0

 

 

0

 

0

 

0

 

y1

y0

hz0

1 0.2 2 1.4

z

2

z

h(2x y

 

z ) 1.6 0.2(2 0.2 1.4 1.6) 1.808;

 

 

 

1

 

1

1

 

1

 

 

y

 

y

hz

1.4 0.2 1.6 1.72

 

2

 

 

 

1

1

 

 

 

 

 

 

 

 

Сведем полученные результаты решения в следующую таблицу: Таблица 6.4-1

xi

yi

zi

0

1

2

0.2

1.4

1.6

0.4

1.72

1.808

Более точное решение ОДУ 2-го порядка можно получить, используя одну из формул Рунге-Кутта 4-го порядка, например,

y

i 1

y

 

h(z

 

h / 6(k

k

 

k

 

));

 

 

i

 

i

1

 

2

 

3

 

zi 1 zi

 

h / 6(k1 2k 22k3 k4 ),

где

k1 f ( xi , yi , zi );

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

k

2

f (x

i

h / 2; y

i

h / 2 z

i

h2

/ 8 k ; z

i

h k

 

/ 2);

 

 

 

 

 

 

 

 

 

 

 

 

 

1

 

 

 

1

 

k

3

f ( x

i

h / 2; y

i

h / 2 z

i

h2

/ 8 k

2

; z

i

h k

2

/ 2);

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

k

4

f (x

i

h; y

i

h z

i

h2

 

/ 2 k

; z

i

h k

3

)

 

 

 

 

 

 

 

 

 

 

 

3

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

56

 

 

 

 

 

 

 

 

 

 

 

 

Для решения ОДУ n-го порядка

y

(n)

 

f(x,y,y ,y ,...,y

(n 1)

)

 

введем следующие

обозначения:

y

(1)

y ', y

(2)

 

 

 

 

В результате этих подстановок порядка:

y '', y(n)

перейдем

y

(n 1)

.

 

к системе из n ОДУ первого

 

 

f (x,y

 

 

,y

 

 

,...,y

 

 

);

y

(1)

(2)

(n 1)

 

1

1

 

 

 

 

 

 

 

y

f

(x,y

 

 

,y

 

 

,...,y

 

 

);

 

(1)

(2)

(n 1)

2

2

 

 

 

 

 

 

 

 

 

 

 

 

 

 

.

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

.

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

y

f

(x,y

(1)

,y

(2)

,...,y

(n 1)

).

 

n

n

 

 

 

 

 

 

 

(6.4-2)

 

 

 

Решением системы (6.4-2)

являются

y

n

φ (x).

 

 

 

 

n

 

 

 

y(n 1)

При

заданных начальных

условиях

y

(n 1)

(x0 ) и использовании метода Эйлера

 

 

 

 

 

 

помощью следующих рекуррентных формул

функции

y

φ (x), y

2

φ (x),...,

 

1

1

2

x x0 , y y(x0 ),..., y(1)

y '(x0 ),...,

решение может быть получено с

y

(i 1)

y

(i)

hf (x

 

i

 

(i)

,y

 

 

(i)

,...,y

 

(i)

);

 

 

,y, y

 

 

 

 

 

 

 

 

 

 

 

 

 

i

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

1

 

1

1

 

 

1

 

 

2

 

 

 

n 1

 

 

y

 

(i 1)

y

(i)

hf (x

 

i

 

(i)

,y

 

(i)

,...,y

(i)

);

 

 

 

 

i

,y, y

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

2

 

 

2

2

 

 

 

1

 

 

 

2

 

 

 

n 1

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

.

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

.

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

y

 

(i 1)

y

(i)

 

 

i

 

(i)

,y

 

(i)

,...,y

(i)

).

n

 

n

hf (xi ,y, y

1

 

2

 

n 1

 

 

 

 

 

n

 

 

 

 

 

 

 

 

 

 

 

Окончательным решением ОДУ n-го порядка, согласно определению, служит функция y φ(x) , вычисленная на заданном множестве точек [a;b].

Тема 7. Одномерная оптимизация

7.1. Постановка задачи

Задача оптимизации – одна из важнейших составляющих многих инженерных задач. Найти оптимальное решение – означает найти наилучший, в смысле заданного критерия, вариант из всех возможных. При решении задачи оптимизации рассматривается некоторая функция, называемая целевой (или

критериальной), и аргументы (параметры целевой функции), называемые параметрами оптимизации.

По количеству независимых переменных различают задачи одномерной оптимизации (n=1) и многомерной оптимизации (n 2). При этом задача нахождения максимума целевой функции сводится к задаче нахождения

57

минимума путем замены функции f(x) на -f(x), поэтому в дальнейшем будем говорить только о поиске минимума функции, то есть такого x* [a, b], при котором f(x*) = min(f(x)).

В области допустимых значений функция f(x) может иметь несколько экстремумов (минимумов или максимумов - рис. 7.1-1). Говорят, что функция f(x) имеет в точке x*локальный минимум, если существует некоторая положительная величина , такая, что если х – х* , то f(x) f(x*), т.е. существует - окрестность точки х*, такая, что для всех значений х в этой окрестности f(x) f(x*). Функция f(x) имеет глобальный минимум в точке x*, если для всех х справедливо неравенство f(x) f(x*). Таким образом, глобальный минимум является наименьшим из локальных минимумов.

y

x1

x

3

x

5

 

 

 

x

 

x4

 

 

2

 

 

 

 

 

 

Рис. 7.1-1

Задачей одномерной оптимизации является нахождение точек локального минимума и соответствующих им значений функции, а в некоторых случаях требуется вычислить глобальный минимум. Однако, во всех случаях эта задача сводится к задаче нахождения локального минимума.

Интервал, на котором локализован единственный минимум, называется

отрезком неопределенности.

Известно, что необходимым условием существования экстремума дифференцируемой функции f(x) является выполнение равенства f (х) = 0. Точка х, удовлетворяющая данному условию, называется точкой стационарности. Достаточным условием существования минимума в точке

стационарности является выполнение неравенства f (х)>0, а максимума - f (х)<0.

Функция является унимодальной на отрезке [a,b], если она на этом отрезке имеет единственную точку глобального минимума и слева от этой точки является строго убывающей, а справа строго возрастающей.

Достаточным условием унимодальности функции f(x) на отрезке [a;b] является следующее: если функция f(x) дважды дифференцируема на отрезке [a,b] и f (х*)>0 в любой точке этого отрезка, то функция f(x) - унимодальна на отрезке [a,b].

58

Заметим, что условие f (х*)>0 определяет множество точек, на котором функция является выпуклой (вниз). Условие f (х*)<0 определяет вогнутую функцию, которая на [a,b] имеет максимум и также является унимодальной.

Все численные методы одномерной оптимизации применяются только для унимодальных функций.

Задача одномерной оптимизации имеет единственное решение в том случае, если функция f(x)на отрезке [a;b] имеет только один экстремум.

Пример 7.1-1.Провести исследование функции f(x) = x3 – x + e-x на предмет существования экстремумов.

Проведем графическое исследование. Построим график функции f(x).

f(x) имеет две точки минимума: х1

и х2, причем точка х1 – точка глобального

минимума. Отрезки неопределенности:

для точки х1 - [-4;-3], а для точки х2 -

[0;1].

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

10

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

5

 

 

 

 

f (x)

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

4

2.5

1

0.5

 

 

 

 

 

 

 

 

 

 

 

5

 

 

 

 

 

 

 

 

 

 

10

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

x

 

 

Проверка достаточного условия существования минимума:

f (x) = 3x2

– 1 – e-x;

f (x) = 6x + e-x,

 

 

f (0)

< 0,

f (1) >0, f (x) > 0 для х 0;1 ,

 

 

f (-4)

< 0,

f (-3) >0, f (x) > 0 для х -4;-3 .

 

 

f (x)

> 0 для всех

х 0;1 и

х -4;-3 . Следовательно, функция f(x)

является унимодальной на выбранных отрезках.

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

Вследующих случаях, когда:

значения функции f(x) определены в ходе эксперимента;

 

целевая функция очень сложна или

не имеет непрерывных

 

производных;

 

 

классические методы поиска оптимального значения не применимы.

 

59

 

Для поиска минимума применяют численные методы одномерной оптимизации, к которым относятся методы одномерного поиска.

Суть методов одномерного поиска заключается в том, что на каждой итерации интервал неопределенности уменьшается и стягивается к точке минимума. Уменьшение отрезка происходит до тех пор, пока на некоторой n-й итерации отрезок неопределенности bn;an не станет соизмеримым с заданной погрешностью , то есть будет выполняться условие |bn-an| < . Тогда за точку минимума можно принять любую точку, принадлежащую этому отрезку, в частности, его середину

X min

a

n

b

 

 

n

 

 

 

2

 

 

 

 

.

 

 

 

 

7.2. Метод прямого перебора с переменным шагом

На практике чаще применяют одну из основных модификаций метода –

метод прямого перебора с переменным шагом. Суть его заключается в следующем. От начальной точки интервала неопределенности [a,b] двигаются с начальным шагом h до тех пор, пока функция в точках разбиения уменьшается (рис.7.2-1). Если функция в очередной точке стала возрастать, то происходит сужение интервала неопределенности путем возврата от рассматриваемой (которая становится правой границей нового интервала) точки x на два шага назад. При этом левой границей нового отрезка неопределенности станет точка a=x-2h, а правой b=x. Новый отрезок вновь исследуют таким же образом, но уже с шагом, уменьшенным в два раза h=h/2. Процесс повторяется до тех пор, пока отрезок неопределенности не станет сопоставимым с заданной точностью

|bn-an| < .

Рис.7.2-1

60

Соседние файлы в предмете Численные методы