ВУЗ: Не указан
Категория: Не указан
Дисциплина: Не указана
Добавлен: 02.08.2019
Просмотров: 4700
Скачиваний: 6
ЛР 7. Метод Эйлера. Схемы Рунге-Кутта решения ОДУ.
141
кет содержит операцию левого матричного деления, с помощью которой
очень просто найти решение: X = A\f. Более того, эта операция позволяет
решать недоопределённые и переопределённые системы линейных уравне-
ний, выбирая алгоритм решения в зависимости от вида матрицы
A.
3
k
Задайте случайным образом матрицу
A размерности 20
× 20 и век-
тор
X. Определите число обусловленности матрицы A с помощью функ-
ции cond. Изменяя значения некоторых элементов матрицы
A, добейтесь,
чтобы её число обусловленности стало больше
10
3
. Используя
A и X, най-
дите вектор
f = AX. Полагая вектор X неизвестным, решите систему
линейных уравнений всеми предложенными выше методами и сравните
найденные решения с уже известным. Какой из методов дал более точный
результат? Обратите внимание на решения, полученные обычным методом
Гаусса и методом с выбором главного элемента.
Реализация функциями Matlab.
\ — операция левого матричного деления
rref(A) — приведение матрицы к диагональному виду
inv(A) — нахождение обратной матрицы
cond(A) — нахождение числа обусловленности матрицы
3.7
ЛР 7. Метод Эйлера. Схемы Рунге-Кутта решения ОДУ.
Рассмотрим обыкновенное дифференциальное уравнение
p-го порядка:
y
(p)
= f (x, y, y
0
, y
00
, . . . , y
(p
−1)
)
Путём введения замены, данное уравнение можно свести к системе линей-
ных уравнений первого порядка:
(
y
k+1
= y
0
k
,
k = 1, 2, . . . , p
− 1,
y
0
p
= f (x, y
1
, y
2
, . . . , y
p
),
ЛР 7. Метод Эйлера. Схемы Рунге-Кутта решения ОДУ.
142
где
y
1
= y. Данная система может быть записана в векторной форме:
dy
dx
= f (x, y)
Для получения единственного решения из системы нужно наложить
p до-
полнительных условий на функции
y
k
(x). Для задачи Коши данные усло-
вия задаются в одной точке:
y
k
(x
0
) = η
k
,
k = 1, 2, . . . , p. Эти условия
рассматриваются как задание начальной точки для интегральной кривой
в
(p + 1)-мерном пространстве (x, y
1
, y
2
, . . . , y
p
). Если правые части систе-
мы непрерывны и ограничены в некоторой окрестности начальной точки
(x
0
, η
1
, η
2
, . . . , η
p
), то решение задачи Коши существует, но может быть не
единственно. Если правые части к тому же удовлетворяют условию Лип-
шица по переменным
y
k
, то решение существует и единственно, т.е. задача
Коши поставлена корректно.
Рассмотрим уравнение 1-го порядка:
(
dy
dx
= f (x, y)
y(x
0
) = η
и пусть данная задача Коши поставлена корректно. Будем искать числен-
ное решение уравнения на отрезке
[x
0
, X]. Введем на этом отрезке сетку
ω
h
=
{x
i
, i = 0, 1, . . . , N
}, таким образом, чтобы x
0
< x
1
< . . . < x
N
= X.
Обозначим
h
i
= x
i+1
− x
i
,
i = 0, 1, . . . , N
− 1 шаг сетки. Заменив произ-
водную в уравнении правой разностью, получим
y
i+1
− y
i
h
i
= f
i
,
i = 0, 1, . . . , N
− 1,
где
f
i
= f (x
i
, y
i
).
Зная
y(x
0
) = η, можно найти все остальные значения y
i
по формуле:
y
i+1
= y
i
+ h
i
f
i
,
i = 0, 1, . . . , N
− 1. Данный метод нахождения числен-
ного решения называется методом Эйлера (или методом ломаных). Схе-
мы, в которых значение функции явно выражается через уже найденные
значения, называются явными, иначе - неявными. Таким образом, схема
Эйлера является явной. Оценка погрешности для данного метода дает
O(max(h
i
)), что предполагает малый шаг сетки для получения удовле-
творительного решения.
ЛР 7. Метод Эйлера. Схемы Рунге-Кутта решения ОДУ.
143
3
k
Найдите численное решение следующего ОДУ методом Эйлера (на
равномерной сетке) и сравните его с аналитическим:
(
dy
dx
= x
2
,
y(0) = 1.
4
k
Matlab имеет множество функций для численного решения обык-
новенных дифференциальных уравнений и их систем. Солверы ode23 и
ode45 основаны на формулах Рунге-Кутты 2,3 и 4,5 порядков соответ-
ственно. Разберем пример их использования на примере задачи о колеба-
ниях под воздействием внешней силы:
(
y
00
+ 2y
0
+ 10y = sin t,
y(0) = 1, y
0
(0) = 0.
Сводим к системе уравнений первого порядка:
y
0
1
= y
2
,
y
0
2
=
−2y
2
− 10y
1
+ sin t,
y
1
(0) = 1, y
2
(0) = 0.
function t e s t ode :
Y0 = [ 1 ; 0 ] ; % вектор начальных условий
[ T Y] = ode45 ( ’ o s c i l ’ , [ 0 1 5 ] , Y0 ) ; % получение решения
% на отрезке 0<t <15
function F= o s c i l ( t , y )
% составляем функцию
F=[y ( 2 ) ;
−2∗y(2) −10∗y(1)+ sin ( t ) ] ;
% из правых частей
вектор Y(:,l) содержит решение исходного уравнения,
вектор Y(:,2) содержит производную решения уравнения.
4
k
Постройте графики координаты
y
1
(t) и скорости y
2
(t). Воспользо-
вавшись знаниями теории обыкновенных дифференциальный уравнений
можно получить аналитическое решение:
y = e
−t
(C
1
cos 3t + C
2
sin 3t) +
1
85
(9 sin t
− 2 cos t),
где для данной задачи Коши
C
1
=
87
85
,
C
2
=
26
85
. Постройте график ана-
литического решения и сравните с численным, полученным при помощи
ode23 и ode45.
ЛР 7. Метод Эйлера. Схемы Рунге-Кутта решения ОДУ.
144
5
k
Решите следующее дифференциальное уравнение:
(
y
00
=
−
1
t
2
,
y(t
0
) = ln t
0
, при t
0
= 0, 01
и сверьте численное решение с аналитическим
y = ln t.
Реализация функциями Matlab. ode45, ode23, ode113, ode15s,
ode23s, ode23t, ode23tb.
Глава 4
Приложения
4.1
Список тем для реферативно-расчётной работы
1. Нахождение всех корней (в том числе комплексных) произвольного
многочлена степени
6
20 методом парабол.
Литература: [4] гл.5, §2, п.8.
2. Интерполяция сплайнами (вычисления методом прогонки)
Литература: [4] гл.2, §1, п.9; [5] гл.3, §4.
3. Интерполяция многочленами Эрмита
Литература: [4] гл.2, §1, п.6; [5] гл.3, §3.
4. Интегрирование методом Симпсона с автоматическим выбором шага
Литература: [5] гл.4, §1, п.5.
5. Решение краевой задачи для дифференциального уравнение 2-го по-
рядка с граничным условием 1-го рода методом прогонки
6. Решение краевой задачи для дифференциального уравнение 2-го по-
рядка с граничным условием 2-го рода методом прогонки
7. Решение краевой задачи для дифференциального уравнение 2-го по-
рядка с граничным условием 3-го рода методом прогонки
8. Решение системы линейных уравнений методом Якоби
9. Решение системы линейных уравнений методом Зейделя
10. Решение системы линейных уравнений методом вращений
145