Введение

I. Интерполирование

§ 1. Интерполирование полиномами

§ 2. Интерполированиие сплайнами

II. Численное интегрирование

§ 1. Квадратурные формулы прямоугольников, трапеций, Симпсона

§ 2. Квадратурные формулы Гаусса

III. Численное решение систем линейных алгебраических уравнений

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

§ 3. Системы с трёхдиогональными матрицами. Метод прогонки

§ 4. Обусловленность систем линейных алгебраических уравнений

§ 5. Итерационные методы решения систем линейных алгебрпических уравнений

IV. Численное решение обыкновенных дифференциальных уравнений

§ 1. Разностные уравнения

§ 2. Численное решение задачи Коши

Литература

    Данное пособие представляет собой систематизированное содержание лекций по курсу “Ведение в численные методы”, который читается на втором курсе факультета ВМиК с 1994 г. и содержит последовательное изложение основных понятий, определений, теорем и утверждений, рассматриваемых и доказываемых на лекциях.
 
 

1.1. Опр. Пусть функция f(х) задана таблично на [a,b]:

x₀ = a, x_n = b , x₀ < x₁ < x₂ < .... < x_n , y_i = f(x_i) i = 0,..., n

Тогда построение непрерывной на [a,b] функции j (x) , такой что j (x_i) = y_i называется интерполяцией функции f(x) на [a,b].

1.2. Опр. Пусть полином степени n L_n(x) = a₀ xⁿ + a₁ x^n-1 + ... + a_n интерполирует y=f(x) на [a,b],
т.е. L_n(x_i) = y_i= f(x_i). Тогда L_n(x) называется интерполяционным полиномом.

Утверждение. Интерполяционный многочлен степени n для функции y=f(x), заданной таблично в n+1 точках, существует и единственен.

Данное утверждение следует из того, что определитель Вандермонда отличен от нуля.
 

1.3. Интерполяционный многочлен в форме Лагранжа имеет вид .

1.4. Теорема.

Пусть функция y=f(x) имеет n+1 непрерывную производную на [a,b], и L_n(x) - интерполяционный многочлен, L_n(x_i) = f(x_i) , i=0,1,...,n. Тогда для погрешности интерполяции y (x) = п L(x) - f(x)п справедлива оценка

  ,

где

1.5. Полиномы Эрмита

Полиномы Эрмита интерполируют таблично заданную функцию с учетом известных значений производной в узлах сетки.

Пусть заданы n+1 узлов x_i , x₀ = a, x_n = b , x₀ < x₁ < x₂ <....< x_n , значения функции в них y_i = f(x_i) и значения производной в них y_i ' = f' (x_i) i = 0,..., n. Требуется построить полином P_2n+1 (x) такой, что P_2n+1(xi) = y_i, P_2n+1' (x_i) = y_i ' . Этот полином и называется полиномом Эрмита.
 
 

Пусть функция y=f(x) задана таблично :

x_{0 =} a, x_{n =} b , x₀ < x₁ < x₂ < .... < x_n , y_i = f(x_i) i = 0,..., n.

2.1. Опр. Кубической сплайн-интерполяцией называется функция j (x) такая, что

    j (x_i) = f(x_i) ,   i=0,1,...,n ,

    j ' (x_i-0) = j ' (x_i+0) ,
                                                                                      (1)
    j ' ' (x_i-0) = j ' ' (x_i+0),    i=1,...,n-1

    j ' ' (x₀) =0, j ' ' (x_n) =0,

и j (x) = a_i + b_i(x-x_i) +c_i(x-x_i)² +d_i(x-x_i)³ , x_i-1 Ј x Ј x_i

    Величины коэффициентов a,b,c,d, находятся из системы уравнений (1). Для нахождения значений этих коэффициентов удобно, с помощью последовательного исключения неизвестных, редуцировать систему (1) к системе трехточечных уравнений относительно коэффициентов с_i , и решать ее далее с помощью метода прогонки (см. далее).
 
 

Опр. Выражение вида 
предназначенное для вычисления определенного интеграла называется квадратурной формулой.

Опр. Величина R=I - S называется погрешностью квадратурной формулы.
 

1.1. Формула прямоугольников

Простая : S₀ = (b-a) f ( (b+a)/2 ) , R₀ = -(b-a)³ f ’’(x)/ 24 , x О (a,b)

Составная: ,  h=(b-a)/n

Формула трапеций

Простая : S₁ = (b-a)( f(b) + f(a) )/ 2, R1 = (b-a)³ f ’’(x)/ 12 , x О (a,b)

Составная: ,  h=(b-a)/n

Формула Симпсона

Простая : S₂ = (b-a)( f (a) + 4f( (b+a)/2 ) + f(b) )/ 6,

R_{2 =} (b-a)⁵ f ⁽⁴⁾(x )/ 90 , x О (a,b)

Cоставная: ,

h=(b-a)/n, здесь - четное число.
 
 

2.1.
В квадратурных формулах Гаусса ищутся не только коэффициенты С_i , но и точки x_i - из соображений обеспечения точности квадратурной формулы для полинома максимальной степени.

Квадратурная формула Гаусса 
будет точна для произвольного полинома степени 2n+1,если величины c_i и x_i удовлетворяют системе уравнений:

2.2.
Можно показать, что узлы x_j квадратурной формулы на отрезке [-1,1] являются корнями полинома Лежандра:

 ,

а коэффициенты квадратурной формулы вычисляются по формулам

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

Найти вектор x , удовлетворяющий уравнению

    Ax = f ,

где A - квадратная матрица порядка n,

    A = { a_{i, j} },  i,j = 1,2,...,n,
    x ^T = {x₁,x₂,..., x_n }, f ^T = { f₁,f₂,...,f_n },

или, что тоже самое, найти x_{1,x2,..., xn} удовлетворяющие системе уравнений

    a ₁₁ x₁ + a ₁₂ x₂ + ... + a _1n x_n = f ₁

    a ₂₁ x₁ + a ₂₂ x₂ + ... + a _2n x_n = f ₂

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

    a _n1 x₁ + a _n2 x₂ + ... + a _nn x_n = f _n
 
 

2.1. Метод Гаусса состоит в приведении матрицы к треугольному виду.

Приведение матрицы к треугольному виду осуществляется по формулам

 ,   m= k+1,...,n, l= k,..., n, k=1,...,n-1.

(прямой ход метода Гаусса), в результате чего получается система

    a⁽ⁿ⁾ ₁₁ x₁ + a⁽ⁿ⁾ ₁₂ x₂ + ... + a⁽ⁿ⁾ _kk x _k + ... + a⁽ⁿ⁾ _1n x_n = j ⁽ⁿ⁾₁

    a⁽ⁿ⁾ ₂₂ x₂ + ... + a⁽ⁿ⁾ _kk x _k + ... + a⁽ⁿ⁾ _2n x_n = j ⁽ⁿ⁾ ₂

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

    a⁽ⁿ⁾ _kk x _k + ... + a⁽ⁿ⁾ _kn x_n = j ⁽ⁿ⁾_k

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

    a⁽ⁿ⁾ _nn x_n = j ⁽ⁿ⁾ _n
 

Вычисление решения x_{1,x2,..., xn} осуществляется следующим образом

(формулы обратного хода Гаусса)

2.2. Общее число действий в методе Гаусса порядка (n^{3 +3n2)/3.}
 
 

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

Метод прогонки применяется для решения систем линейных уравнений с матрицей специального вида - трехдиагональной матрицей :

    - a_i x_i-1 + c_i x_i - b_i x_i+1= f_i,    i=2, ... , n-1

    x₁ = k ₁ x₂ + n ₁ (2)

    x_n = k ₂ x_n-1 + n ₂
 

На первом этапе находятся коэффициенты

a _i+1 = b_i /(c_i - a_i a _i) , b _i+1 = ( b _i a_i + f_i ) /(c_i - a_i a _i ) i=2, ... , n-1 (прямой ход), a ₂ = k ₁ , b ₂ = n _1,

а на втором этапе находится решение

    x_i = a _i+1 x_i+1 + b _i+1 , i = n-1, ... ,2,1

    x_n = (n ₂ + b _n ) / ( 1-a _n k 2 ).
 

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

Теорема

Пусть система уравнений (2) такова, что   a_i > 0, b_i > 0, c_i > 0, c_i і a_i + b_i , i=2,..., n-1,
и пk ₁п + пk ₂п <2, пk ₁п Ј 1, пk ₂п Ј 1.

Тогда метод прогонки корректен.

Теорема

Пусть система уравнений (2) такова, что a_i > 0, b_i > 0, c_i > 0, c_i і a_i + b_i , i=2,..., n-1,
и существует i₀ >1 такое, что c_i0 > a_i0 + b_i0 , и п k ₁п Ј 1, п k ₂п Ј 1.

Тогда метод прогонки корректен.

Теорема

Пусть система уравнений (2) такова, что зa_iз , зb_i з ,зc_i з >0, зc_i з > зa_iз+...+зb_i з , i=2,..., n-1,
и п k ₁п Ј 1, п k ₂п Ј 1.

Тогда метод прогонки корректен.
 
 

4.1. Опр. Под нормой матрицы А будем понимать 
 

5.1. Канонический вид одношаговых итерационных методов для решения системы линейных алгебраических уравнений
Ax = f :

где В_к+1 - квадратная матрица nґ n , t _k+1 > 0 - итерационный параметр. В дальнейшем будем использовать следующие согласованные нормы в конечномерном пространстве размерности n:

- евклидова норма 

- норма в С 

- энергетическая норма A=A^*>0

5.2. Итерационный метод сходится, если  .

Опр. Величина z_k = x_k - x называется погрешностью решения.

Опр. Если B_k+1 = B и t _k+1 = t то метод называется стационарным
 

Теорема.

Пусть A=A^* >0, и B- 0.5t A>0 тогда итерационный метод

сходится в норме кк * кк _A , .
 

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

Канонический вид метода простой итерации - 

Теорема

Пусть A=A^*>0, , тогда метод простой итерации сходится.
 

Замечание. Если A=A^*>0 то у матрицы А существует n собственных значений l_k таких, что

    0< l _min=l _1<l 2<...<l _n-1<l _n=l _max ,

при этом и .

Можно найти оптимальное значение итерационного параметра t = t ₀ такое, что заданная точность решения будет достигаться за минимальное число итераций.

Теорема
Пусть A=A^*>0, тогда при для погрешности явного метода простой итерации z_k справедлива оценка .
 

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

Каноническому виду метода Зейделя соответствует B=(D+L), t =1, где D -диагональная матрица, L - нижняя треугольная матрица 
 

Индексный вид метода Зейделя: 

Теорема

Пусть A=A^*>0, тогда метод Зейделя сходится.

Теорема

Пусть матрица A такова, что, i=1,..., n, к qк <1, тогда метод Зейделя сходится со скоростью геометрической прогрессии со знаменателем q,
т.е. .
 

Метод релаксации

Канонический вид:

w < 1- метод нижней релаксации

w = 1- метод Зейделя

w > 1- метод верхней релаксации

Теорема Пусть A=A^*>0, 0<w <2, тогда метод релаксации сходится.
 
 

1.1. Опр. Пусть дан отрезок [a,b]. Равномерной сеткой на этом отрезке назовем множество узлов w _h такое, что  w _h = { x_j = jh,  j=0,...,n, h=(b-a)/n) }.

Опр. Сеточной функцией y = y_j = y(x_j) называется функция, заданная в узлах сетки.

Любую сеточную функцию y _j = y(x_j) можно представить в виде вектора
Y=(y₀, y₁, ..., y_n-1, y_n), и, следовательно, множество сеточных функций образует конечномерное пространство, в данном случае размерности n+1. В этом пространстве можно ввести норму, например
  или .

Пусть дано дифференциальное уравнение
Lu(x) = f(x,u) ( например, ) .

1.2. Заменим Lu в узле сетки x_i линейной комбинацией значений сеточной функции y_i на некотором множестве узлов сетки, называемом шаблоном. Такая замена Lu на L_h y_h называется аппроксимацией на сетке дифференциального оператора L разностным оператором L_h . Замена непрерывной функции f(x,u) в узлах сетки на сеточную функцию j (x_h,y_h) называется аппроксимацией правой части.

Таким образом дифференциальное уравнение можно аппроксимировать (заменить) на сетке разностной схемой

L_h y_h = j ( x_h,y_h) ( например, ).

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

Пусть u_h - проекция непрерывной функции u(x) на сетку ( например, u_h = u(x_j) = u_j ).

Опр. Погрешностью аппроксимации дифференциального оператора Lu разностным оператором L_h назовем величину
y ₁ = (Lu)_h - L_h u_h , где (Lu)_h - проекция на сетку результата действия дифференциального оператора L на функцию u

( например, ) .

Опр. Говорят, что погрешность аппроксимации дифференциального оператора имеет в узле x_i порядок k , если
y ₁(x_i) = O(h^k) ® 0 при h® 0.

Опр. Погрешностью аппроксимации правой части f сеточной функцией j _h назовем величину y ₂ = f_h - j _h , где f_h - проекция на сетку функции f(x,u) (например, f(x_{j ,uj).}

Опр. Погрешность аппроксимации правой части имеет в узле x_i порядок m , если y ₂ = O(h^m) ® 0 при h® 0

Опр. Погрешностью аппроксимации разностной схемы на решении в узле x_i (локальной погрешностью) назовем величину y , равную

y = y ₁ - y ₂ = (Lu)_h - L_h u_h - ( f_h - j _h )= j _h - L_h u_h ,

здесь использовано, что Lu=f.

.

Опр. Значения локальной погрешности аппроксимации в каждом узле x_i образуют сеточную функцию погрешности аппроксимации y _i .

Обычно требуется оценка погрешности аппроксимации на сетке, т.е. оценка функции y _i в некоторой сеточной норме.

Опр. Говорят, что погрешность аппроксимации разностной схемы
имеет m-ый порядок на сетке, если чкyчк = O(h^m).

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

ч к z_h ч к =ч к u_h - y_h ч к = O(h^k) ® 0 при h® 0.
 
 

2.1. Требуется найти численное решение задачи Коши для ОДУ первого порядка

На отрезке [0,T] вводим разностную сетку w t = { t_n = nt , n=0,...,N, t =T/N } и сеточную
функцию y_n = y(t_n).
 

_{2.2. Метод Эйлера}

Явный метод Эйлера.

Погрешность аппроксимации явного метода Эйлера y = O(t ) - первого порядка
( y ₁ = O(t ), y ₂ = 0 ).

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

Величины p_k и a _k выбираются из соображений аппроксимации.
 

Метод Рунге - Кутта второго порядка

Величины p₁, p_{2 ,} a определяются из условия второго порядка аппроксимации,
для данной схемы p₁ = (2 a -1)/ a , p₂ = 1/2 a .

Метод Рунге- Кутта второго порядка можно записать в другом виде

 
 
Условие второго порядка аппроксимации (1-s )a =0.5.
 

На практике широко распространен Метод Рунге - Кутта четвертого порядка :
 
 

Данная схема имеет четвертый порядок аппроксимации.
 
 

2.4.
Другую возможность для повышения порядка сходимости разностных схем для решения ОДУ предоставляют Методы Адамса.

В методах Адамса коэффициенты b_i находятся из условий наивысшего для данного p порядка погрешности аппроксимации.

Методы Адамса являются многошаговыми и являются разновидностью p - шаговых методов. P-шаговый метод может быть записан в следующем общем виде:

P- шаговому методу Адамса соответствует набор коэффициентов a _i вида

a₀ =1, a₁ = -1, a₂ = ... = a_p=0
 
 

Явный двухшаговый метод Адамса.
 

Погрешность аппроксимации имеет второй порядок. Для нахождения y₁ (со вторым порядком) используется обычно схема Рунге-Кутта второго порядка.
 
 

Неявный двухшаговый метод Адамса.
 

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

[1] А.А. Самарский. Введение в численные методы. Москва:  Наука, 1987.
[2] А.А. Самарский, А.В. Гулин. Численные методы. Москва:  Наука, 1989.
[3] А.Н.Тихонов, Д.П.Костомаров. Вводные лекции по прикладной математике. Москва: Наука, 1984.