ВУЗ: Не указан
Категория: Не указан
Дисциплина: Не указана
Добавлен: 02.08.2019
Просмотров: 4724
Скачиваний: 6
Численное решение дифференциальных уравнений
96
Коэффициенты
a
i
,
b
ij
,
σ
i
выбираются из соображений точности. Например, для то-
го, чтобы уравнение (11.10) аппроксимировало исходное уравнение (11.2), необходимо
потребовать
m
P
i=1
σ
i
= 1. Отметим, что методы Рунге–Кутта при m > 5 не используются.
Для получения очередного значения
y
k+1
в методе Рунге–Кутта тре-
буется выполнить много промежуточных вычислений. Но из-за высокого
порядка погрешности аппроксимации метода можно выбирать сетку с до-
вольно крупным шагом
h. В итоге общее число арифметических операций
уменьшается. Наиболее распространённым является метод Рунге–Кутта
четвёртого порядка:
y
k+1
= y
k
+
h
6
(p
1
+ 2p
2
+ 2p
3
+ p
4
),
p
1
= f (x
k
, y
k
),
p
2
= f (x
k
+
h
2
, y
k
+
h
2
p
1
),
p
3
= f (x
k
+
h
2
, y
k
+
h
2
p
2
),
p
4
= f (x
k
+ h, y
k
+ hp
3
).
Глава 2
Практические занятия
2.1
Занятие 1
2.1.1
Ошибки в вычислениях, числа с плавающей точкой
1.1
Пусть надо решить систему двух линейных уравнений:
−10
−7
x
1
+ x
2
= 1,
x
1
+ 2x
2
= 4.
Решение. Первый метод. Исключая
x
1
из первого уравнения:
x
1
= 10
7
x
2
−
10
7
, и подставляя это выражение во второе уравнение, получаем
x
2
=
10
7
+4
10
2
+2
.
Проведя вычисления с семью значащими цифрами, получаем
x
2
= 1.000000,
x
1
= 0.000000, что совершенно неверно, как видно из второго уравнения.
Второй метод. Исключая
x
1
из второго уравнения:
x
1
= 4
− 2x
2
,
получаем для
x
2
формулу
x
2
=
1+4
·10
−7
1+2
·10
−7
. После вычислений получаем
x
2
= 1.000000, x
1
= 2.000000 — правильное (с точностью по шести де-
сятичных цифр) решение.
Представление чисел с плавающей точкой в Matlab По умолча-
нию переменные имеют тип «плавающая точка двойной точности». Каж-
дое число занимает 8 байт. Представление имеет следующую структуру:
8 байт
z
}|
{
b
63
|{z}
знак
s
b
62
. . . b
52
|
{z
}
порядок
e
b
51
. . . b
0
|
{z
}
мантисса
f
97
Занятие 1
98
число с пл. точкой
=
(
−1)
s
(2
e
−1023
)(1.f ) нормализованное, 0 < e < 2047,
(
−1)
s
(2
e
−1022
)(0.f ) ненормализованное, e = 0, f > 0,
ошибка
иначе
.
Далеко не все действительные числа точно представимы числом с пла-
вающей точкой. Множество
F чисел с плавающей точкой имеет мощность
|F |
6
2
64
, в то время как
|
R
| = ∞. Расстояние между числом x ∈ F и бли-
жайшим к нему числом
y
∈ F (x < y) приблизительно пропорционально
|x|. Чем дальше от нуля, тем реже встречаются числа из F среди чисел
из
R
.
Функция Matlabeps(x) вычисляет расстояние до ближайшего числа
с плавающей точкой, большего x.
1.2
Зная внутренне представление чисел системы Matlab, определить,
чему равно eps(7) и eps(8). Ответы сравнить и объяснить, почему они
отличаются в два раза.
Решение. Имеем
7
10
= 111
2
= 1, 11
2
· 2
2
= (
−1)
0
(2
1025
−1023
)(1, 11). Отсю-
да
s = 0, e = 1025, f = 11 00 . . . 0
|
{z
}
50 нулей
2
. Следующее число
7
0
с плавающей
точкой, большее 7, будет иметь
s = 0, e = 1025, f = 11 00 . . . 0
|
{z
}
49 нулей
1
2
. Полу-
чаем, что
eps(7) = 7
0
− 7 = 2
−50
≈ 8.881784197001252 · 10
−16
. Аналогично
eps(8) = 8
0
− 8 = 2
−49
≈ 1.776356839400251 · 10
−15
.
Ответ eps(7) в два раза меньше eps(8). Так происходит всякий раз
когда переходим через пороговое значение вида
2
k
,
k
∈
Z
. В нашем случае
7 < 2
3
, а
8 = 2
3
.
1.3
Объяснить, почему при вычислении в среде Matlab
100 + 2^
−50
получается ответ 100, но
1 + 2^
−50
Занятие 1
99
возвращает ответ 1.000000000000001?
Решение. Найдем двоичное представление первой суммы:
100 + 2
−50
=
1100100
2
+ 0, 0 . . . 0
| {z }
49 нулей
1
2
= 1100100, 0 . . . 0
| {z }
49 нулей
1
2
. Получилось 57 бит. Ман-
тисса имеет длину 52 бита, то есть данное число необходимо округлить
(отрезать «лишние» 5 бит). После округления остаётся только первое сла-
гаемое 100.
Найдем двоичное представление второй суммы:
1+2
−50
= 1
2
+0, 0 . . . 0
| {z }
49 нулей
1
2
=
1, 0 . . . 0
| {z }
49 нулей
1
2
. Получилось 51 бит, что вполне умещается внутри мантис-
сы.
Замечание. 52 двоичных разряда соответствуют 15–16 десятичным.
2.1.2
Матричные вычисления в Matlab
Полезные функции:
• ones(m,n) и zeros(m,n) — создать матрицу размера m × n, все эле-
менты которой 1 и 0. Пример:
ones(2,2)
=
1 1
1 1
!
,
zeros(1,3)
=
0 0 0
.
• eye(n) — создать единичную матрицу размера n × n. Пример:
eye(3)
=
1 0 0
0 1 0
0 0 1
.
• diag(v,k) — расположить элементы вектора v на k-ой диагонали
квадратной матрицы соответствующего размера. Пример:
v
= [1 2]
diag(v,0)
=
1 0
0 2
!
,
diag(v,1)
=
0 1 0
0 0 2
0 0 0
,
diag(v,-1)
=
0 0 0
1 0 0
0 2 0
.
Занятие 1
100
Пусть далее
A =
1 2 3
4 5 6
!
.
• flipud(A) (от англ. flip up down), fliplr(A) (от англ. flip left right),
’ (штрих) — перевернуть матрицу сверху вниз, слева направо, относи-
тельно главной диагонали. Пример:
flipud(A)
=
4 5 6
1 2 3
!
,
fliplr(A)
=
3 2 1
6 5 4
!
,
A’
=
1 4
2 5
3 6
.
• prod(A,k), sum(x) — перемножить, просуммировать элементы матри-
цы в направлении
k. Пример:
prod(A,1)
=
4 10 18
,
prod(A,2)
=
6
120
!
,
sum(A,1)
=
5 7 9
,
sum(A,2)
=
6
15
!
.
• triu(A,k) (от англ. upper triangle), tril(A,k) (от англ. lower triangle)
— вернуть верхнюю, нижнюю треугольную часть матрицы
A, начиная
с
k-ой диагонали. Пример:
triu(A,1)
=
0 2 3
0 0 6
!
,
tril(A,0)
=
1 0 0
4 5 0
!
.
• repmat(B,m,n) — продублировать матрицу B по вертикали m раз и
по горизонтали
n раз. Пример:
B =
1 4
repmat(B,2,3)
=
1 4 1 4 1 4
1 4 1 4 1 4
!
.
1.4
Пусть задан вектор x = [1,2,4,5,6,7]. Задать в Matlabматрицу
определителя Вандермонда:
∆ =
1 x
0
x
2
0
· · · x
n
0
1 x
1
x
2
1
· · · x
n
1
... ...
. . . ...
1 x
n
x
2
n
· · · x
n
n
.