ВУЗ: Не указан
Категория: Не указан
Дисциплина: Не указана
Добавлен: 02.08.2019
Просмотров: 4719
Скачиваний: 6
Занятие 7
126
Для этого построим в явном виде обратную матрицу.
A
−1
=
1 1 2 4
· · · 2
n
−3
2
n
−2
0 1 1 2
· · · 2
n
−4
2
n
−3
... ... ... ... ...
...
...
0 0 0 0
· · ·
1
1
0 0 0 0
· · ·
0
1
,
||A
−1
||
∞
= 1 + 1 + 2 + 2
2
+ . . . + 2
n
−2
= 2
n
−1
,
||A||
∞
= n и µ
∞
(A) = n2
n
−1
,
т.е. матрица
A плохо обусловлена, хотя det(A) = 1.
2.7
Занятие 7
2.7.1
Решение дифференциальных уравнений
7.1
Проверить, аппроксимирует ли разностная схема уравнение
y
0
(x) = f (x, y(x))
а)
1
3h
(y
k
− y
k
−3
) = f
k
−1
;
б)
1
8h
(y
k
− 3y
k
−2
+ 2y
k
−3
) =
1
2
(f
k
−1
+ f
k
−2
);
в)
1
2h
(3y
k
− 4y
k
−1
+ y
k
−2
) = f
k
.
7.2
Для задачи
y
0
+ y = x + 1, y(0) = 0 рассматривается схема
y
k+1
− y
k
−1
2h
+ y
k
= kh + 1,
y
0
= 0,
y
1
= 0.
Каков порядок аппроксимации у данной схемы? Можно ли его улучшить?
7.3
Для задачи
y
0
+ 5y = 5, y(0) = 2 построена разностная схема
y
k+1
− y
k
−1
2h
+ 5y
k
= 5,
y
0
= 2,
y
1
= 2
− 5h.
Исследовать её аппроксимацию.
Занятие 7
127
7.4
Построить аппроксимацию второго порядка для заданной краевой
задачи
u
00
(x) + u(x) = cos x + 1,
u(1) = 2,
u
0
(3)
− 3u(3) = 1,
Решение. Как видно из граничных условий, решение ищется на отрезке
[1, 3]. Разобьём этот отрезок на n равных частей. Длина каждого частич-
ного отрезка равна
h = (3
− 1)/n. Введём в рассмотрение сетку ω =
{x
0
, x
1
, . . . , x
n
}, где x
0
= 1, x
1
= 1+h, x
2
= 2+2h, . . . , x
n
−1
= 3
−h, x
n
= 3.
Рассмотрим также сеточную функцию
y
i
= y(x
i
), определённую только в
узловых точках
x
i
,
i = 0, 1, . . . , n.
Требуется построить систему линейных уравнений с неизвестными
y
i
,
где
y
i
≈ u(x
i
). Совокупность таких y
i
будет численным решением краевой
задачи.
Погрешность аппроксимации должна по условию задачи иметь порядок
2 (т.е., например, уменьшение шага h вдвое должно уменьшать погреш-
ность решения
|y
i
− u(x
i
)
| в четыре раза).
Заменим в уравнении производную разделённой разностью (см. лек-
цию 5)
u
00
(x
i
)
≈
y
i
−1
−2y
i
+y
i+1
h
2
. Исходное дифференциальное уравнение пе-
рейдёт в разностное
u
00
(x) + u(x) = cos x + 1
→
y
i
−1
− 2y
i
+ y
i+1
h
2
+ y
i
= cos x
i
+ 1.
Убедимся, что получен второй порядок аппроксимации. В разностное урав-
нение вместо сеточной функции подставим точное решение, т.е. заменим
y
i
на
u
i
= u(x
i
). Кроме этого разложим u
i
±1
в ряд Тейлора в окрестности
точки
x
i
u
i
±1
= u
i
± hu
0
i
+
h
2
2
u
00
i
±
h
3
3!
u
000
i
+
h
4
4!
u
IV
(ξ
±
), ξ
−
∈ [x
i
−1
, x
i
], ξ
+
∈ [x
i
, x
i+1
].
После сокращения имеем
u
00
i
−
h
2
12
[u
IV
(ξ
−
) + u
IV
(ξ
+
)] + u
i
= cos x
i
+ 1.
Занятие 7
128
Вычитая из последнего равенства исходное дифференциальное уравнение,
получаем невязку
ψ
(1)
n
=
−
h
2
12
[u
IV
(ξ
−
) + u
IV
(ξ
+
)] =
−
h
2
6
u
IV
(η) = O(h
2
), где
η
∈ [ξ
−
, ξ
+
]. Получен требуемый порядок аппроксимации.
Очевидно, замена граничного условия
u(1) = 2 на y
0
= 2 не имеет
погрешности аппроксимации.
В правом граничном условии нельзя заменить
u
0
(3) левой разделённой
разностью, так как
u
0
(3) =
y
n
−y
n
−1
h
+ O(h
1
). Выделим в O(h
1
) главный
член. Из формулы Тейлора в окрестности
x
n
= 3 имеем
u(3
− h) = u(3) − hu
0
(3) +
h
2
2
u
00
(3) + O(h
3
),
откуда
u
0
(3) =
u(3)
− u(3 − h)
h
+
h
2
u
00
(3) + O(h
2
).
Из исходного уравнения следует, что
u
00
(3) + u(3) = cos(3) + 1. Таким
образом,
u(3)
− u(3 − h)
h
+
h
2
[cos(3) + 1
− u(3)] + O(h
2
)
− 3u(3) = 1.
Невязка для правого граничного условия равна
O(h
2
), что даёт второй
порядок аппроксимации.
Окончательный ответ — это система из
n + 1 линейного уравнения с
неизвестными
y
0
, y
1
, . . . , y
n
y
i
−1
− 2y
i
+ y
i+1
h
2
+ y
i
= cos x
i
+ 1, где i = 1, 2, . . . , n
− 1 и x
i
= 1 + ih,
y
0
= 2,
y
n
− y
n
−1
h
+
h
2
[cos(3) + 1
− y
n
]
− 3y
n
= 1.
Глава 3
Лабораторный практикум
3.1
ЛР 1. Распространение ошибок в вычислительных проце-
дурах.
При построении математической модели и получении результата, на ошиб-
ку полученного ответа влияют несколько факторов:
1. Ошибка математической модели — ошибка неучтенных параметров
физического явления. Например, полет тела, брошенного под углом к
горизонту с малыми скоростями можно описывать обычными балли-
стическими уравнениями без учёта сопротивления воздуха. Как толь-
ко скорости становятся значительными нужно вводить сопротивление
воздуха. Далее следуют: учёт плотности, температуры воздуха, вет-
ра, переменного ускорения свободного падения, вращения Земли и т.д.
Построение адекватной математической модели описывается в курсе
математического моделирования.
2. Ошибка входных данных — это ошибка измерений параметров физи-
ческой модели. Во многих физических и технических задачах данная
погрешность достигает нескольких процентов. Она приводит к так на-
зываемой неустранимой погрешности, которая не управляется мате-
матически и, зачастую, приводит к парадоксальным результатам.
3. Ошибка численного метода — связана с тем, что исходные операторы
заменяются приближёнными. Например, интеграл – суммой, диффе-
ренцирование – конечной разностью, функцию – полиномом, беско-
нечный ряд – конечной суммой элементов. Погрешность численного
129
ЛР 1. Распространение ошибок в вычислительных процедурах.
130
метода обычно берут в 2-5 раза меньше неустранимой погрешности.
Меньше брать невыгодно из-за увеличения объёма выполняемых вы-
числений, больше — из-за снижения точности вычислений. Ошибка
численного метода управляется математически.
4. Ошибка округления связана с ограниченной разрядной сеткой ком-
пьютера, из-за чего, например, бесконечные десятичные дроби пред-
ставляются в конечном виде.
3
k
Начнём со второго пункта. Рассмотрим многочлен Уилкинсона
P (x) = (x
− 1)(x − 2)...(x − 20) = x
20
− 210x
19
+ 20615x
18
+ . . .
Очевидно, корнями его являются
x
1
= 1, x
2
= 2, . . ., x
20
= 20. Выполните
в Matlab команду p=poly(1:20), которая позволяет получить коэффи-
циенты полинома, корнями которого является аргумент функции. Набе-
рите roots(p) и убедитесь, что корни полинома найдены верно. Измени-
те значение, например, второго коэффициента на малую величину
10
−7
(составляет
≈ 5 · 10
−8
% от числа) и снова выполните эту команду - по-
ловина корней стала комплексными числами! Исходная задача оказалась
неустойчивой к входным данным (их малое изменение ведет к сильному
изменению решения), в результате чего получился абсурдный результат.
3
k
Рассмотрим представление числа в компьютере. Чаще всего в Mat-
lab производятся операции над числами с двойной точностью. Число
формата double — 64-разрядное число (8 байт), в котором 1 бит — знако-
вый, 52 отводятся под мантиссу и 11 под порядок числа. (
D =
±(1+m)·2
n
,
m = 0, m
1
m
2
. . . m
k
,
m
1
6= 0 — мантисса числа, n — порядок). Так как в
порядке тоже есть знаковый бит, то множество его значений лежит от
−2
10
+ 1 =
−1023 до 2
10
= 1024. В командном окне Matlab набери-
те 2ˆ1023. Получилось число с десятичным порядком
+308 (его легко
оценить аналитически:
2
1024
≈ 2
1000
= 1024
100
≈ (10
3
)
100
= 10
300
). Те-
перь выполните 2ˆ1024 — получили Inf т.е. машинную бесконечность.
Команда realmax как раз и выдаёт максимальное число, которое можно
представить в Matlab, команда realmin — минимальное (по абсолютной