Введение. Численные методы и приближённые вычисления

Лекция 1

 

Введение. Численные методы и приближённые вычисления

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

Главная задача прикладной математики нахождение решений с требующейся точностью. Этим она отличается от классической математики, которая основное внимание уделяет исследованию условий существования решений и его свойств.

В истории прикладной математики можно различить 3 периода:

1 период 3-4 тыс. лет назад – вычисление объемов, площадей и решение простых задач арифметики, алгебры, геометрии. (жрецы)

2 период начался с Ньютона. Решались задачи астрономии, геодезии и расчета механических конструкций. Задачи сводились либо к решению ОДУ, либо систем линейных алгебраических уравнений (СЛАУ). Бруно, Коперник.

3 период начался с 1940 года, когда надо было быстро рассчитывать, быстро передвигаться, предметы (катюша), зенитные орудия и т.д.

Курчатов, Сахаров, Королев, Появились ЭВМ.

Это потребовало от исследователей разработки и совершенствования ранее существующих методов и разработки новых методов.

Не следует думать, что знание численных методов и новых ЭВМ позволяет сразу решить любую задачу.

Решение сложной инженерной задачи выполняется в следующей последовательности.

Формулируется конечная цель решения, выявляются исходные и получаемые параметры.

Система отношений, связывающая исходные и конечные параметры.

Позволяет по исходным данным найти исходный результат. Составляется схема последовательных действий (алгоритм).

 

Тот же алгоритм, но записанный на языке «понятном» для ЭВМ.

 

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

Раздел вычислительной математики состоит из следующих частей:

1) Общие понятия об интерполировании функции. Численное интегрирование и дифференцирование.

2) Решение СЛАУ.

3) Решение нелинейных уравнений.

4) ОДУ.

5) Задачи на собственные решения.

Методы решения

Методы решения задач делятся на аналитические, графические и численные.

Численные методы- это методы приближенного или точного решения математической задачи, основанные на построении конечной последовательности арифметических действий над числами. Численные методы связанны, прежде всего, с погрешностями при работе с числовой информацией.

Абсолютная и относительная погрешность.

Пусть α- точное значение некоторой величины, α_n- известное приближенное к этому точному значению. Абсолютной погрешностью приближения α называется величина

или α_n= α ± ∆α

Абсолютная погрешность далеко не всегда может охарактеризовать эффективность метода получения приближённого результата: например

∆α= 1 мм при измерении расстояния в 1 км или диаметра вала 4 мм, имеют совершенно различную характеристику качества измерения. Поэтому часто пользуются относительной погрешностью

 

или

Часто используют относительную погрешность выраженную в %.

Лекция 2

Линейная интерполяция.

При линейной интерполяции предполагается, что функция f(x) между узлами интерполяции изменяется по линейному закону (см. рис)

 

из геометрии уравнение прямой, проходящей через точки (x_i-1 y_i-1) и (x_i y_i) будет иметь вид (

Отсюда для каждого x, лежащего в интервале [x_i-1 x_i ] может быть найдено соответствующее значение y по формуле

Таким образом, мы можем найти любое f(x) =y на любом отрезке интерполяции.

Лекция 3

Методы приближенного интегрирования.

Формулы Ньютона –Котеса.

Один из методов приближенного интегрирования f(x)d(x) заключается в замене функции f(x) на какой- либо интерполяционный полином. Например

В этом случае, если взять узлы интерполяции равноотстоящими, то

x_k=x₀+k*h k=0,1,2…n

x=x₀+t*h dx=hdt

x-x₀=(x₀+th)-(x₀+0h)=th

x-x₁=(x₀+th)-(x₀+1h)=(t-1)h

_ _ _ _ _ _ a=x₀=x₀+th t=0

x-x_n=(x₀+th)-(x₀+nh)=(t-n)h

 

x_k-x₀=(x₀+kh)-(x₀+0h)=kh

_ _ _ _ _ b=x_n=x₀+nh=x₀+nh t=n

x_k-x_k-1=(x₀+kh)-(x₀+(k-1)h)=1h

x_k-x_k+1=(x₀+kh)-(x₀+(k+1)h)=-1h

_ _ _ _ _ _

x_k-x_n=(x₀+kh)-(x₀+nh)=-(n-k)h

 

dt* f(x_k)= A_n^(k)f(x_k)h

… dx= (b-a)ⁿ

[a, b]

Частные случаи.

n=1 L₁(x)=ax+ b

 

A₁⁽⁰⁾=- dt=

A₁⁽¹⁾= dt=

f(x)dx=h()=

f(x)dx= = h[ ]

n=2 L₂(x)=a² x+ b+c

A₂⁽⁰⁾=- dt=

A₂⁽¹⁾= dt=

A₂⁽²⁾= dt=

 

f(x)dx=h()=

f(x)dx= = = [ ]

 

Лекция 4

 

Лекция 5

 

 

Разделённые разности

 

Пусть функция f(x) вычислена в точках i=0,1,2,…,n, . Тогда выражение

называется разделённой разностью 1-го порядка.

Разность вида: называется разделённой разностью 2-го порядка

- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

Разность вида: -называется разделённой разностью к-го порядка. Вычисляются разделённые разности согласно Таблицы №1.

 

 

Таблица №1

 

1) Разделённые разности симметричны относительно своих аргументов

т. е. узлы можно менять любым образом.

2) Разделённые разности порядка k от являются однородным многочленом относительно своих аргументов степени n-k, при n=k она равна 1, а при k>n она равна 0.

 

1

0

1 0

0

1

 

 

3) Линейность разделённой разности относительно функции, к которой применяется разделённая разность.

 

 

Конечные разности.

В случаях когда узлы интерполяции берутся равноотстоящими

x_k=x₀+k*h k=0,1,2…n

n- шаг интерполирования, то часто вводят следующие обозначения

– конечные разности

 

0нисходящие

восходящие - конечная разность 2 ^го порядка

центральные

- конечная разность 3 ^го порядка

 

Лекция 6

 

Отделение корней.

Пусть [a, b] – интервал изменения х, в котором отыскиваются вещественные корни уравнения f(x)=0, где f(x) – непрерывная функция вещественного переменного. Разбив [a, b ] на n отрезков длинной , наличие простого корня в интервале [x_k-1, x_k], k=1, 2, …, n, легко установить по знаку произведения концевых значений функции f:

Алгоритм определения корней можно показать следующей блок схемой (см. лист 9). Алгоритм предусматривает ввод (блок 1) исходных данных (a, b – начало и конец исследуемого интервала, n – число одинаковых отрезков разбития [a, b] на подынтервалы) и определение величины (блок 2). Далее для каждого интервала длиной вычисляется (блок 3) и y=f(x_k) (блок 7), а затем, в результате исследования знака произведения (блок 8), решается вопрос о существовании на простого корня. При предусматривается выдача значений границ (блок 9) и выявление начала следующего интервала (блок 10).

 

Если отказаться от блока 10, то при случайном совпадении корня с точкой x_k деление отрезка произойдёт выдача на печать значений двух интервалов и , хотя нам достаточно одного.

При определении корней алгебраического уравнения (*) можно воспользоваться следующими закономерностями:

а) алгебраические уравнения п – го порядка имеет п корней, среди которых могут быть вещественные и комплексные;

b) число положительных вещественных корней равно числу (или меньше на чётное число) перемен законов в последовательности коэффициентов a₀, a₁ , a₂, …, a_n, причём равные нулю коэффициенты не учитываются (теорема Декарта).

 

 

После отделения корней можно приступать к нахождению корня с заданной точностью на каждом выделенном интервале.

Метод итераций.

Приведём уравнение к виду (1)

Возьмём некоторую произвольную точку на оси . Если существует такая точка , что , то очевидно, число и будет корнем уравнения (1). В общем случае . Суть метода простых итераций (повторений) состоит в построении последовательности { } точек, сходящихся к корню , где , (2)

При этом называют начальным приближением.

 

Необходимо отметить, что последовательность { } не всегда является сходящейся. Рассмотрим пример:

Линейное уравнение , имеющее корень можно преобразовать к виду или .

Взяв за начальное приближение , в первом случае получим

; , …

Данная последовательность стремится к точному решению уравнения.

Во втором случае

; ; , …

Эта последовательность расходится.

Следовательно, вид функции оказывает существенное влияние на сходимость итерационного процесса.

Сходимость метода простой итерации определяется следующей теоремой: При любом начальном приближении на некотором интервале последовательность , построенная в соответствии с методом простой итерации (2), сходится к корню , который является единственным решением уравнения (1) на рассматриваемом интервале, если здесь функция удовлетворяет условию Липшица

с постоянной Липшица .

Если функция дифференцируема и имеет на ограниченную производную, то постоянная Липшица .

В рассмотренных выше примерах для уравнения и итерационный процесс является сходящимся к корням уравнения .

Для уравнения и условия сходимости не выполняется.

Если вычисление постоянной Липшица сопряжено с определёнными математическими трудностями, то алгоритм метода простой итерации использует косвенную оценку точности вычислений

,

где - допустимая погрешность вычислений.

Блок – схема алгоритма метода простой итерации с косвенной оценкой точности вычислений:

 

 

Лекция 7

 

Численные решения СЛАУ.

Метод Гаусса.

Пусть задана система линейных уравнений Ax=U с плотной матрицей, элементы которой хранятся в оперативной памяти.

 

 

Переходим к системе (исключения х₁)

Эта операция равносильна умножению Ax=U слева на матрицу

 

где p_i1= i=2,3…n A₁ x= U₁à Ax= U

a_ij⁽¹⁾= a_ij- a_ij

Следующий переход к системе (исключения х₂)

Это исключения равносильно умножению A₁ x= U ₁ слева на матрицу Р₂:

где p_i2= i=3,4…n a_ij⁽²⁾= a_ij⁽¹⁾ - a_2j⁽¹⁾

A₁ x= U₁=> A₁х = U₁àP₂P₁Ax= P₂P₁ U

Продолжая этот процесс и так далее, мы придем к легко решаемой системе уравнений с треугольной матрицей A_n-1х =U_n-1

 

 

A_n-1 =

 

A_n-1- верхняя треугольная матрица. Здесь A_n-1х=U_n-1=> =P_n-1P_n-2Ax= P_n-1P_n-2 U, так что A =(P_n-1P_n-2 …P₁)^-1 A_n-1х

 

P Q

 

A = *

 

Где матрица P- нижняя треугольная, а Q- верхняя треугольная. Решение по методу Гаусса сводится к разложению A на произведение двух треугольных матриц. Одна из них

P~ с единицами по главной диагонали и на верхнюю треугольную

 

Q~

 

Заметим, что

det A=det P*det Q= a₁₁ a₂₂⁽¹⁾ a₃₃⁽²⁾… a_nn^(n-1)

Если представления матрицы A в виде PQ найдены, то решение системы Ax=U сводится к решению двух треугольных систем уравнений.

Py= U ↓ (y)

Ax=U=>PQx=U=>

Qx= y ↑ (x)

Решения по методу Гаусса сводится фактически к решению системы Py= U ↓

При прямом ходе, а затем к решению системы Qx= y ↑ при обратном ходе.

Разложение матрицы на две треугольные матрицы называется процессом факторизации (триангуляция) матрицы, т.е. представление матрицы в виде A=PQ.

После вычисления прямого хода мы получаем матрицу A_n-1= Q

Метод Холецкого, прогонки приводят к факторизации матрицы A.

Метод Холецкого.

Пусть дана невырожденная матрица размером N*N (det A≠0) Представим её в виде произведения

A=BC (1) A~(a_ij) B~(b_ij) C~(c_ij)

 

 

k>i

b_ik=0, когда k>i

c_kj=0 при k>j c_jj =1

Из (1) следует

a_ij= b_ik c_ki

Преобразуем эту сумму двумя способами

a_ij= b_ik c_ki= b_ik c_ki+ b_ii c_ij+ b_ik c_ki= b_ik c_ki+ b_ii c_ij

 

a_ij= b_ik c_kj= b_ik c_kj+ b_ij c_jj+ b_ik c_kj= b_ik c_kj+ b_ij c_jjà b_ij c_jj= a_ij- b_ik c_kj

Отсюда находим

b_ij= i ≥ j b₁₁=a₁₁ c_jj=1

c_jj= при I < j

Матрицы B и C найдены, тогда Ax= BCx=U

By=U ↓ (y) прямой ход

 

Cx=y ↑ (x) обратный ход

 

Лекция 8

Метод прогонки.

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

Пусть требуется найти решение следующей системы трехточенных уравнений.

 

Или в векторной форме

A (2)

A=

По методу Гаусса A= * = *

Факторизация ленточных матриц, т.е. разложение их на произведение двух треугольных матриц A=PQ

(нижней и верхней) треугольных матриц может быть произведено так, что P и Q также будут ленточными матрицами.

Перейдем к следующей системе (3)

y_j=p_j+1y_j+1+q_j+1 (3)

Где p_j+1 и q_j+1 неизвестные числа

u₀- p₁u₁ =q₁

u₁- p₂u₂ =q₂ =

u₀- p₁u₁ =q₃

 

 

Числа p_j и q_j можно искать из системы (1), т.к.

y_j-1=p_jy_j+q_j (4)

Подставляя (4) в систему (1) получим

-a_j(p_j y_j+ q_j)+c_j y_j-b_j y_j+1=f_j

(c_j-a_j p_j)y_j=b_j y_j+1 +(f_j+ a_j q_j)

(5)

Сравнивая (5) с(3) получим

(6)

Из 1 ^го крайнего условия

y₀ =p₁y₁+q₁

=> (7)

 

Из 2 ^го крайнего условия

y_n-1= p_ny_n+q_n

-a_n y_n-1+c_ny_n=f_{n (8)}

Теперь с помощью (3) получим

Проверка метода прогонки на устойчивость.

Поскольку вычисления p_j и q_j производятся по рекуррентным формулам, то в зависимости от элементов матрицы A: a_i, b_i, c_i, то неизбежны ошибки округления при вычислении по рекуррентным формулам p_j и q_j. Эти ошибки могут быстро возрастать и сделать вычисления непригодными. Поведение ошибки dp_j+1 можно проследить, положив её приближённо равной дифференциалу функции.

Имеем

(9)

Прогонка будет устойчивой, если и (10)

Тогда dp_j+1 <dp_j

Эти неравенства (10) будут выполнены, если, например, положить, что коэффициент матрицы A удовлетворяют следующие условия

0≤ a_j≤b_j

c_j≥ a_j+b_j (11) j=

0<

В самом деле

, отсюда

<

Вообще, проверяется p_j+1≤ 1 для j=1,n-1

При этих условиях при вычислении p_j+1 накопление ошибки не происходит. Аналогично из равенств

;

<dq_j

Т.к. при тех же условиях (11) нарастания погрешностей при вычислении q _j+1 не происходит.

 

Лекция 9

Лекция 10

 

 

Сплайны.

Пусть есть функция f(x) заданная в узлах отрезка [a,b]. Пусть заданны узлы интерполяции x₀x₁…x_n x_i [a,b] в которых заданны значения функции f(x₀) f(x₁)… f(x_n) [x₀x₁], [x₁x₂]… [x_n-1x_n]~∆- разбиения.

Будем приближать функцию f(x) на каждом отрезке ∆~[x_j-1x_j] полиномом

P_j,m(x)=a_0j+a_1jx+…+ a_njx^m

Сплайном назовем S_∆(x) функцию совпадающую на каждом участке [x_j-1x_j] с полиномом P_j,m(x) и такую, что во всех внутренних узлах непрерывную вместе со своими производными до (m-1) порядка включительно во всех точках x_i i=1, n-1.

Кубический сплайн.

Будем считать, что на каждом участке [x_j-1x_j] мы приближаем полином 3=ей степени.

P_j,3(x)=ax³+bx²+cx+d, тогда

(1)

Найдем S_∆(x) для чего проинтегрируем дважды (1), получим

Найдем с₁ и с₂ из условия S_∆(x_j-1)= y _j-1 S_∆(x_j)= y _j

C₁h_j=y_j- M_j => C₁h_j=y_j-y_j-1+

C₁=

C₂=

Поскольку S_∆^`(x_j-0)= S_∆(x_j+0)

Приравнивая получаем

 

(2)

λ_j= μ_i= λ_j+ μ_i=1

μ_iM_j-1+2M_j+ λ_j M_j=d_j j=1,n-1

Полученная система служит для определения M₀ M₁… M_n, подставляя в формулы (х) и (х^*) получим, что полином удовлетворяет 3-м свойствам

1) S_∆(x_j-0)= S(x_j+0)

2) S<sup>`</sup>(x<sub>j</sub>-0)= S` (x_j+0)

3) S^{``</sup>(x<sub>j</sub>-0)= S<sup>``}(x_j+0)

 

Лекция 11

Метод простой итерации.

Пусть задана система линейных уравнений

Ax=b=>Ax-b=0 (1) Ax`-b=ξ- невязка

x- x+ Ax-b=0=>x= (I-A) x+ b=>x=Bx+ b (2)

Пусть матрица B невырожденная, т.е. det B≠0/

Выбрав произвольно x₀ вычисляем последовательно x_k

x_k= Bx_k-1+b (3)

Пусть x^*- точное решение уравнения (1) и (2), тогда

 

(x_k-x^*) =B (x_k-1-x^*) (4)

(x_k-x^*) =B (x_k-1-x^*) = B B (x_k-2-x^*) = B^k (x₀-x^*) (4a)

Пронормируем последнее выражение (4а) получим:

(5)

Из полученного неравенства следует, что при произвольном выборе x₀ соотношение (6) при < 1 (7)

Эквивалентно или . Соотношение (6) возможно, если выполнено неравенство (7). < 1

Лекция 12

Метод Зейделя.

Решение системы Ax=U по методу Зейделя производится по формулам

 

Если при вычислении i-той координаты вектора

учитывается найденные заранее уже координаты …

То вычисления будут проходить по формулам Зейделя.