Численные методы решения инженерных задач

МИНИСТЕРСТВО ОБРАЗОВАНИЯ И НАУКИ РОССИЙСКОЙ ФЕДЕРАЦИИ ФЕДЕРАЛЬНОЕ ГОСУДАРСТВЕННОЕ БЮДЖЕТНОЕ ОБРАЗОВАТЕЛЬНОЕ УЧРЕЖДЕНИЕ ВЫСШЕГО ОБРАЗОВАНИЯ «ТУЛЬСКИЙ ГОСУДАРСТВЕННЫЙ УНИВЕРСИТЕТ» Институт горного дела и строительства Кафедра ССМиК Теличко Виктор Григорьевич доцент, кандидат технических наук КОНСПЕКТ ЛЕКЦИЙ по дисциплине ЧИСЛЕННЫЕ МЕТОДЫ РЕШЕНИЯ ИНЖЕНЕРНЫХ ЗАДАЧ «Численные методы решения инженерных задач» Уровень профессионального образования: высшее образование – бакалавриат Направление подготовки: 08.03.01 «Строительство» Профиль подготовки: «Промышленное и гражданское строительство» «Городское строительство и хозяйство» Квалификация выпускника: бакалавр Форма обучения: заочная Тула 2015 1 Рассмотрено на заседании кафедры протокол № Зав. кафедрой от « » 20 г. А.А.Трещев 2 СОДЕРЖАНИЕ 1. Лекция 1 2. Лекция 2 3. Лекция 3 4. Лекция 4 5. Лекция 5 6. Лекция 6 7. Лекция 7 8. Лекция 8 9. Лекция 9 3 ВВЕДЕНИЕ Лекция 1 ВВЕДЕНИЕ В ТЕОРИЮ ЧИСЛЕННЫХ МЕТОДОВ ЦЕЛЬ ЛЕКЦИИ: определить предмет курса; дать общую характеристику таких свойств численных методов как эффективность, точность, устойчивость, экономичность; пояснить эти свойства на простейших вычислительных правилах. Основные требования, предъявляемые к вычислительным алгоритмам Вычислительный алгоритм должен быть устойчивым к ошибкам, обеспечивать требуемую точность, быть эффективным по времени выполнения, быть экономичным по требуемым объемам памяти, не допускать аварийных остановов. Поясним эти требования. Устойчивость. Пусть последовательность величин yi , i  1,2 , вычисляется по рекурсивному правилу yi 1  yi  d, i  0 ,1,2 , при заданных y0 и d. Предположим, что при вычислении yi внесена ошибка i (например, за счет операции округления), т. е. вместо значения yi использу-ется приближенное значение ~ yi  yi   i . В соответствии с рекурсивным правилом ~ y i 1  ~ y i  d  y i   i  d  y i 1   i . Следовательно, ~ yi 1  yi 1   i , и ошибка, допущенная на i-м шаге процесса, на следующем шаге не увеличивается, если операция сложения выполнена без новых округлений. Это означает, что алгоритм устойчив. Рассмотрим другое рекурсивное правило вычисления этой последовательности: yi 1  qyi , i  0,1,2, Опять y0 и q считаем заданными. Пусть, как и в предыдущем примере, ~ y i  yi   i . Тогда ~ y i 1  q~y i  q( y i   i )  qy i  q i , или ~ yi 1  yi 1  q i . В этом случае ~ yi 1  yi 1  q i , и при q>1 ошибка будет возрастать. Такой алгоритм является неустойчивым. Точность. Обозначим через y точное решение задачи, yk – ее приближенное значение, полученное на k-м шаге численного метода. Тогда y  y  yk составляет погрешность решения, а   y  y  y k является абсолютной погрешностью. Вычислительный алгоритм должен давать решение с заданной точностью  1 и, следовательно, критерием завершения процесса уточнения решения является выполнение неравенства   y  y k  1 . 4 На практике в силу трудностей вычисления абсолютной погрешности вместо  используют ее оценку сверху – предельную абсолютную погрешность Δy. В качестве Δy выбирают как можно меньшее значение, удовлетворяющее неравенству y   , а критерием завершения процесса в этом случае является неравенство  y  1 . Часто используют относительную погрешность   y, т. е. абсолютную погрешность на единицу измерения, и предельную относительную погрешность y  . Критерий завершения процесса в этом случае имеет вид  y  2 , где  2 – заданная допустимая относительная ошибка. Рассмотрим в качестве примера, как может быть построена оценка предельной абсолютной погрешности вычисления значений функции u  f  x1,x2 , ,xn  , если известны абсолютные погрешности xi ее аргументов. Ошибка имеет вид u  f  x1  x1,x2  x2 , ,xn  xn   f x1,x2 , ,xn  . Тогда n f n f Δ  u   xi   xi , i 1 xi i 1 xi и в качестве предельной относительной погрешности можно использовать величину n  f Δxi i 1 xi . y   y1 ,y2 , ,yn  Δu  Пусть решением задачи является вектор  T , который будем рассматривать как элемент векторного пространства ние Rn . Приближенное реше- T ~ y  ~ y1,~y2 , ,~ yn  и его погрешность  y  ~ y  y1  ~ y1,y2  ~ y2 , ,yn  ~ yn  T также являются элементами этого пространства. Рассмотрим, как оцениваются ошибки решения в этом случае. Для этого используем количественную характеристику вектора в виде нормы. Говорят, что в R n задана норма, если каждому вектору y из R n сопоставляется вещественное число y , называемое нормой вектора y , для которого справедливы следующие свойства: 1. y  0 , причем y  0 тогда и только тогда, когда y  0 , 2. y   y для любого вектора y и любого числа , 3. y  z  y  z для любых векторов y и z . В вычислительных методах наиболее употребительны следующие нормы: 12 n  n 2 y   yi ; y    yi  ; y  max yi 1 2  1 i  n i 1 i 1  . Абсолютную и относительную погрешность вектора в любой из перечисленных выше норм можно определить следующим образом: ~ y   y  ~ y ;  ~ y  y y  ~ . y 5 Имеет смысл говорить об абсолютной и относительной погрешности (mn)- матрицы решений A. В этом случае используется такая характеристика, как норма матрицы, согласованная с нормой вектора. Для нормы матрицы A, обозначаемой A , справедливы следующие свойства: 1. A  0 , причем A  0 тогда и только тогда, когда A  0 , 2. A   A для любой матрицы A и любого числа  , 3. A  B  A  B для любых (mn)-матриц A и B . Каждой из векторных норм соответствует своя согласованная норма матриц. В частности, нормам y 1 , y 2 , y  соответствуют нормы A 1, A 2 , A  , вычисляемые по правилам m  A 1  max  a ij , A 2  max  j A A 1 j n i  1 1 j n T , n A   max  a ij , 1im j  1 где j(ATA) – собственные числа матрицы A A. Абсолютная и относительная погрешности матрицы решений вычисляются через соответствующие нормы: ~ A  A ~ ~  ~  A  A A ;  A  A   , ~ причем A  – искомая матрица решений, A – ее некоторое приближение. Эффективность. Эффективность вычислительного алгоритма оценивают по времени реализации либо по общему числу операций, требуемых для получения результата. Обозначим через nсл число операций сложения, через nу – число операций умножения, через nдел – число операций деления, через tу, tдел, tсл – время выполнения операций умножения, деления и сложения на ЭВМ. Объективный критерий сравнения вычислительных алгоритмов – требуемое время вычислений: T=nуtу+nслtсл+nделtдел. Если tу≈tдел≈tсл, то T≈(nу+nдел+nсл)tоп, где tоп – время выполнения на ЭВМ арифметической операции, и критерием оценки эффективности алгоритма может быть nΣ=nу+nдел+nсл . Если tу≈tдел>>tсл, то критерий эффективности – число так называемых длинных операций nдл=nдел+nу. Пусть Q – критерий эффективности метода. Значение Q зависит от требуемой точности . Из двух сравниваемых вычислительных алгоритмов решения задачи из них будет эффективнее тот, который при заданном значении ε характеризуется меньшим значением Q. В качестве примера построим два алгоритма вычисления полинома n-й степени Pn  x   an x n  an 1x n 1    a1x  a0 . Непосредственная реализация: P0  a0 , P1  a1 x  a0  P0  a1 x, P2  a 2 x 2  a1 x  a0  P1  a 2 x 2 ,  Pn  Pn 1  a n x n . Расчет полинома по рекурсивному правилу Pk  Pk 1  ak x k ,k  1,2, ,n , вытекающему из приведенной реализации, соответствует алгоритму 1. Алгоритм 1. 1. Начальная установка P0  a0 , e0  1. 2. Вычислить в цикле (k=1,2,…,n): ek  ek 1x, Pk  Pk 1  ak ek . 6 Алгоритм 1 характеризуется таким числом арифметических операций: nу=2n; nсл=n. Преобразуем исходную запись полинома: Pn x   (( (an x  an 1 ) x  an  2 ) x    a1 ) x  a0 . Теперь расчет полинома можно выполнить следующим образом: bn  an , bn 1  bn x  an 1, bn  2  bn 1 x  an  2 ,  b0  b1 x  a0  Pn  x  . Записанное правило реализуется алгоритмом 2. Алгоритм 2. 1. Начальная установка: bn  a n . 2. Вычислить в цикле (k=n,n-1,…,1): bk 1  bk x  a k 1 . Значение b0  Pn x  . Затраты на выполнение алгоритма: nу=n , nсл=n. Очевидно, что алгоритм 2 эффективнее по времени выполнения алгоритма 1. Экономичность. Экономичность алгоритма оценивается по требуемым для его реализации объемам памяти ЭВМ. Аварийные остановы. Любое число ЭВМ принадлежит не всей числовой оси, а некоторому интервалу M 0 ,M   . Если какое-либо число x в процессе работы алгоритма выходит за указанный интервал, то происходит так называемый аварийный останов. Необходимо построить вычислительный алгоритм таким образом, чтобы вычисленные в процессе его работы числа x  M 0 ,M   . Часть 2. СИСТЕМЫ ЛИНЕЙНЫХ АЛГЕБРАИЧЕСКИХ УРАВНЕНИЙ Лекция 2 ПРЯМЫЕ МЕТОДЫ ЦЕЛЬ ЛЕКЦИИ: Определить два класса численных методов (прямые и итерационные); показать, как строятся прямые методы Гаусса, LU-факторизации, Холесского; выполнить оценку их эффективности. Постановка задачи. Основная задача вычислительной алгебры – решение систем линейных алгебраических уравнений (СЛАУ) a11 x1  a12 x 2    a1n x n  b1, a 21 x1  a 22 x 2    a 2n x n  b2 ,  a n1 x1  a n 2 x2    a nn x n  bn . В дальнейшем будем использовать запись этой системы в компактной форме: n  a ij x j  b i ; i  1,n j 1 ( запись i  1,n означает, что индекс i изменяется от 1 до n с шагом 1), или в векторном виде Ax  b , где  a11 a A   21     a n1 a 12 a 22  a n2     a 1n   b1   x1      a 2 n  b x2 2 ; b    ; x   .              a nn   b n   x n  7 Предполагается, что матрица A неособенная, т. е. det A  0 , и решение единственно. Прямые и итерационные методы. Численные методы решения СЛАУ делятся на две большие группы: прямые и итерационные. Прямые методы при отсутствии ошибок округления за конечное число арифметических операций позволяют получить точное решение x  . В итерационных методах задается начальное приближение x 0 и строится последовательность x k x , k   где k – номер итерации. В действительности итерационный процесс прекращается, как только x k становится достаточно близким к x  . Имеется промежуточный класс методов, в которых решение ищется итерационно, однако для них заранее известно, какое число итераций необходимо выполнить, чтобы в отсутствии ошибок округления получить точное решение. На практике при вычислении приближенного решения число итераций в наиболее эффективных методах оказывается значительно меньшим, чем этого требует теория точного решения. Какой класс методов лучше? Однозначно на этот вопрос ответить нельзя. Итерационные методы привлекательнее с точки зрения объема вычислений и требуемой памяти, когда решаются системы с матрицами высокой размерности. При небольших порядках системы используют прямые методы либо прямые методы в сочетании с итерационными методами. Метод Гаусса. В методе Гаусса линейная система Ax  b решается в два этапа. На первом этапе система Ax  b преобразуется к виду (см. рис. 2.1) Ux  q , Рис. 2.1. Структура системы и портрет ее ненулевых элементов до (а) и после (б) прямого хода Гаусса где U – верхняя треугольная матрица с единичной диагональю (это так называемый прямой ход Гаусса). На втором этапе (обратный ход Гаусса) решается система Ux  q . Рассмотрим эти этапы подробнее. Прямой ход. Прямой ход Гаусса состоит из n шагов. Первый шаг. Полагаем, что a11  0 и разделим на него первое уравнение. Перепишем систему с учетом этого преобразования: a1 j  b , j  2,n ; q1  1 ,  x1  u12 x2    u1n xn  q1; u1 j  a a  11 11  a x    a x  b , i  2 ,n . in n i  i1 1 Умножим первое уравнение на ai1 и вычтем его из i-го уравнения преобразованной системы: (ai 2  ai1u12 ) x2    (ain  ai1u1n ) xn  bi  ai1q1, i  2 ,n. Обозначим aij(1)  aij  ai1u1 j , bi(1)  bi  ai1q1 . Получим x1  u12 x2    u1n xn  q1; (1) ai(21) x2    ain xn  bi(1) , i  2 ,n. Второй шаг. На втором шаге из системы (1) ai(21) x2    ain xn  bi(1) , i  2,n исключается x2 аналогичным образом: 8 x2  u23 x3    u2 n xn  q2 ; ( 2) ai(32) x3    ain xn  bi( 2) ; i  3,n, (1) (1) u 2 j  a2(1j) a22 ; j  3,n; q2  b2(1) a22 ; aij( 2)  aij(1)  ai(12)u 2 j ; bi( 2)  bi(1)  ai(21) q2 ; i  3,n; j  3,n. K-й шаг. Запишем общий вид преобразованной системы после k-го шага прямого хода Гаусса: x1  u12 x2  u13 x3    u1k xk  u1, k 1xk 1    u1n xn  q1, x2  u 23 x3    u 2k xk  u2, k 1 xk 1    u2 n xn  q2 ,  xk  uk , k 1 xk 1    ukn xn  qk , (k ) ai(,kk)1 xk 1    ain xn  bi( k ) , i  k  1,n. Здесь ( k 1) ( k 1) ( k 1) ( k 1) ukj  akj akk ; j  k  1,n; qk  bk(k 1) akk ; akk  0, ( k 1) ( k 1) aij(k )  aij(k 1)  aik u kj ; bi( k )  bi( k 1)  aik qk , i  k  1,n; j  k  1,n. Проиллюстрируем, как меняется матрица системы в процессе прямого хода Гаусса на примере системы четвертого порядка (рис. 2.2; ненулевые элементы матрицы обозначены крестиками). Рис. 2.2. Преобразование матрицы системы 4-го порядка на прямом ходе Гаусса Оценим количество длинных операций (умножений и делений) на первом шаге прямого хода Гаусса. Преобразование первого уравнения требует n таких операций. Преобразование остальных n-1 уравнений – n(n-1) операций умножения и деления. Таким образом, первый шаг выполняется за n(n  1)  n  n 2 длинных операций. Рассуждая по аналогии, нетрудно найти затраты на остальных n-1 шагах. Суммарные затраты прямого хода Гаусса определяются в итоге рядом пр. 2 2 2 2 nдл .  n  ( n  1)    2  1  n(n  1)(2n  1) n 3  6 3 . Последняя оценка имеет место для n>>1. Обратный ход. Запишем систему, решаемую на обратном ходе, в координатном виде x1  u12 x2  u13 x3    u1,n 1xn 1  u1n xn  q1, x2  u 23 x3    u2 ,n 1 xn 1  u 2n xn  q2 ,  xn 1  u n 1,n xn  qn , xn  q n . Ее решение: n xn  qn ; xk  qk   u kj x j ; k  n  1,1. j  k 1 Запись k  n  1,1 означает, что индекс k изменяется от значения n-1 до 1 с шагом 1. Требуемое число длинных операций на обратном ходе обр. nдл .  1  2  3    (n  1)  n(n  1) n 2  . 2 2 9 Приближенная оценка справедлива для n>>1. Общие затраты метода Гаусса: пр. обр. n  nдл .  nдл.  nn  12n  1 nn  1 n3   . 6 2 3 Таким образом, при больших n основные затраты в методе Гаусса приходятся на прямой ход. Метод LU-факторизации. В методе LU-факторизации (эту схему называют компактной схемой Гаусса) при решении системы Ax  b выполняется следующая последовательность действий. Матрица A представляется в виде произведения A  LU , Рис. 2.3. Структура матриц L и U в разложениях Дулиттла (а) и Краута (б) где L  нижняя треугольная матрица, U  верхняя треугольная матрица. Такое разложение единственно при условии предварительного выбора диагональных элементов одной из матриц. В этом случае число элементов в матрице A совпадает с суммарным числом неизвестных элементов матриц L и U. Если диагональ L принимается единичной, то такое разложение называют разложением Дулиттла (рис. 2.3,а), если единична диагональ U – разложением Краута (рис. 2.3,б). В дальнейшем при построении метода LU-факторизации будем привлекать разложение Краута. Система Ax  b заменяется системой LUx  b , легко решаемой за два шага: Шаг 1. Ly  b . Принимая во внимание треугольный вид матрицы L, нетрудно получить, что в алгоритме Краута y1  Шаг 2. k 1  1  b1 ; yk  b k   l kj y j ; k  2,n.  l 11 l kk  j 1  Ux  y . Решение этой системы в алгоритме Краута: n  u kj x j , k  n  1,1. . x n  y n , xk  y k  j  k 1 Суммарные затраты реализации обоих шагов при n>>1 составляют n 2 длинных операций. Получим соотношения для расчета элементов матриц L и U в алгоритме Краута. Для этого перемножим матрицы L и U и приравняем результат к A. По правилу перемножения матриц n aij   lik u kj . k 1 Учтем, что lik  0 , k  i; u kj  0 , k  j; u kj  1 , k  j. Рассмотрим элемент aij (рис. 2.4), расположенный на центральной диагонали либо в нижней треугольной части матрицы A. Для такого элемента i ≥ j. Из рисунка следует, что 10 Рис. 2.4. Иллюстрация вычисления элемента матрицы, расположенного ниже главной диагонали j j 1 k 1 k 1  lik u kj lij   lik u kj , aij  так как i ≥ j и u jj  1 . Отсюда j 1 lij  aij   lik u kj , i  j. k 1 Рассмотрим элемент aij (рис. 2.5), находящийся выше главной Рис. 2.5. Иллюстрация вычисления элемента матрицы, расположенного выше главной диагонали диагонали матрицы A (для него j>i). В этом случае i i 1  lik ukj  lii uij   lik u kj . aij  k 1 k 1 Следовательно, uij  1 lii i 1    aij   lik ukj  , i  j.   k 1   Получили в итоге соотношения, которые позволяют вычислять элементы матриц L и U. Последовательность вычислений: сначала вычисляется столбец матрицы L, далее строка матрицы U, затем опять столбец матрицы L, далее строка матрицы U и т. д. (см. рис. 2.6, который иллюстрирует последовательность вычислений и схему хранения матриц L и U). Вычисление столбца матрицы L и строки матрицы U назовем шагом LU-разложения. Приведем в качестве примера схему хранения элементов матриц A,L,U после второго шага LUразложения (рис. 2.7). Число длинных арифметических операций на этапе LU-разложения при n>>1 составляет величину n 3 3 , на шаге решения линейных систем с треугольными матрицами – n 2 . Суммарное число длинных операций приближенно равно n 3 3 (как и в методе Гаусса), а) a11 a21 a31 a41 б) a12 a22 a32 a42 a13 a23 a33 a43 a14 a24 a34 a44 l11 l 21 l 31 l 41 в) u12 l 22 l 32 l 42 u13 u 23 l 33 l 43 u14 u 24 u 34 l 44 1 5 6 2 8 11 12 3 9 13 15 4 10 14 16 7 Рис. 2.6. Исходная матрица A (а), схема хранения L и U матриц (б), последовательность вычисления элементов в принятой схеме хранения (в) 11 l11 l 21 l 31 l 41 u12 l 22 l 32 l 42 u13 u 23 a33 a43 u14 u24 a34 a44 Рис. 2.7 Схема хранения элементов 44 матриц A, L, U после второго шага LU-факторизации т. е. основные затраты приходятся на LU-факторизацию матрицы A. Эта особенность делает особо привлекательным метод LU-факторизации при решении СЛАУ с одной и той же матрицей A, но разными правыми частями: Ax  b i , i  1,2, ,m ; m  n. В этом случае факторизация матрицы выполняется однократно, требуя n3 3 длинных операций, а решение каждой системы с соответствующей правой частью реализуется за n 2 таких операций. Метод Холесского. Исключительно эффективную реализацию метода LU-факторизации можно получить, если ограничиться классом линейных систем с симметрической положительно определенной матрицей A, т. е. A  A T  0. Такую реализацию называют методом Холесского, либо методом квадратного корня. Будем полагать, что решаемая система Ax  b имеет симметрическую положительно определенную матрицу A. В этом случае матрица A представляется в виде ~~ A  LL T . ~ Здесь L – нижняя треугольная матрица. Такое разложение существует и единственно для положительно определенных симметрических матриц. Система преобразуется к виду ~~ LL T x  b . Вектор x ищется путем последовательного решения двух систем с треугольными матрицами: ~ ~ Ly  b ; L Т x  y . ~ Для получения расчетных соотношений элементов матрицы L рассмотрим произвольный элемент матрицы A: j ~ ~ aij   lik l jk , i  j. k 1 Суммирование здесь выполняется только до j, т. к. j≤i. Выделим член при значении k=j: j 1~ ~ ~~ aij  lij l jj   lik l jk . k 1 Теперь ~ 1 lij  ~ l jj j 1~ ~   aij   lik l jk  , i  j; k 1   j 1~ ~ l jj  a jj   l jk2 , i  j. k 1 ~ Эти соотношения позволяют вычислить по столбцам элементы матрицы L . Эффективность такого метода достигается на этапе разложения матрицы, т. к. необходимо ~ вычислить в этом случае только матрицу L . Арифметические затраты в методе Холесского составляют n 3 6 длинных операций и n операций извлечения квадратного корня. Существует другой вариант разложения симметрической положительно определенной матрицы, в котором удается избежать операций извлечения квадратного корня. В этом варианте вводится новая матрица L по правилу 12 ~ ~ L  L diag l jj  , причем ~ L ~ – матрица, вычисленная ранее по схеме Холесского, diagl jj  – диагональная матрица с ~ ~ элементами l jj матрицы L . Матрица L существует и является нижней треугольной с единичной диагональю. В этом случае ~ ~ ~~ A  LLТ  L diag l jj diag l jj LТ  L D LТ , где ~ D  diagl jj2 . Расчетные соотношения для элементов матриц L и D можно получить, как и прежде, привлекая правило перемножения матриц j aij   l ik d kk l jk , i  j , k 1 из которого следует, что j 1 lij d jj  aij   lik d kk l jk , i  j, k 1 j 1 d jj  a jj   lik d kk lik , i  j , k 1 т. к. матрица L имеет единичную диагональ. Такой алгоритм потребует вдвое большего числа перемножений, чем схема Холесского. Однако, если ввести замену переменных aij  lij d jj , то расчетные соотношения примут вид j 1 aij  aij   aik l jk , k 1 j 1 d jj  a jj  i  j,  aik lik , i  j. k 1 Здесь сначала вычисляют вспомогательные величины aij , а затем их используют для определения искомых величин d jj и lij . Количество умножений при такой организации алгоритма составляет приблизительно n 3 / 6 и не содержит операции извлечения квадратного корня. Лекция 3 УСТОЙЧИВОСТЬ И ТОЧНОСТЬ ПРЯМЫХ МЕТОДОВ ЦЕЛЬ ЛЕКЦИИ: Выполнить оценку устойчивости и точности прямых методов; показать, как перестановкой строк и столбцов обеспечивается устойчивость и точность прямых методов, каким образом осуществляется выбор матриц перестановки строк и столбцов в случае систем с разреженной матрицей. Устойчивость алгебраических методов. Прямые методы в отсутствии ошибок округления позволяют получить точное решение системы Ax  b . Современные вычислительные машины оперируют с конечными десятичными дробями, представленными в форме с плавающей точкой. В этом случае уже на этапе запоминания элементов матрицы A и вектора b в памяти ЭВМ вносится ошибка округления и реально решается возмущенная система 13 A  F x  b  h (1) (1) . Здесь F и h – возмущения матрицы системы и вектора правой части. Для элементов f ij(1) матрицы F (1) и компонент hi(1) вектора h (1) справедливы оценки (1) (1) 1 f ij    aij , hi1   bi , где aij – элемент матрицы A, bi – компонент вектора b , ε – число, характеризующее относительную погрешность машинной ариф-метики. Например, для двоичных ЭВМ, использующих арифметику с плавающей точкой и t – разрядную мантиссу,   2 t . Перейдем к более общей числовой оценке возмущений – норме. Из записанных выше неравенств следует, что (1) 1 F  A , h  b , где знаком  обозначена какая-либо норма вектора и согласованная с ней норма матрицы. Поясним теперь суть одного из наиболее разработанных подходов к анализу устойчивости алгебраических методов. Пусть ~x – приближенное решение СЛАУ, полученное применением некоторого прямого метода. Очевидно, что вследствие ошибок округления при реализации на ЭВМ прямого метода это решение не будет точно удовлетворять системе Ax  b . Пусть, однако, удается показать, что ~x является точным решением системы A  F ~x  b  h . Если для матрицы F и вектора h , называемых эквивалентными возмущениями метода, можно получить оценки вида F  f n  A , h  hn  b , где f(n), h(n)  некоторые степенные функции типа cn k с небольшим показателем k, то метод считается устойчивым по Уилкинсону. Такой вид функций f(n), h(n) объясняется тем, что в процессе реализации прямых методов на ЭВМ неизбежно некоторое накопление ошибок округления, пропорциональное числу арифметических операций, за которое прямой метод приводит к решению. Точность алгебраических методов. Рассмотрим влияние возмущений матрицы системы и вектора правой части на точность решения. Относительная погрешность решения  x  ~x  x  x  , где x  – точное решение, в любой согласованной норме подчиняется неравенству ~ x x  x   A A 1  F h .    b  1  A 1 F  A Обозначим   A A 1 ,  A  F A , b  h b . Величины δA , δb являются оценками относительных возмущений матрицы A и вектора b. Число χ есть не что иное, как число обусловленности матрицы A. Естественно допустить, что A 1 F  1 . Такое неравенство означает, что норма матрицы возмущений существенно меньше нормы исходной матрицы, т. е. предполагается, что прямой метод решения линейной системы устойчив. В этом случае  x    A   b  . Относительная ошибка решения системы прямым методом может достигать величины, определяемой произведением числа обусловленности и суммы относительных возмущений 14 матрицы системы и вектора правой части. Если число обусловленности велико, то даже при малых эквивалентных возмущениях матрицы и вектора правой части погрешность решения линейной системы оказывается значительной. Задачи, для которых χ велико, называют плохо обусловленными. Число обусловленности зависит от выбора матричной нормы. Например, для симметричной положительно определенной матрицы A в качестве A можно взять ее максимальное собственное число. Тогда число обусловленности   max min , где max , min – максимальное и минимальное собственные числа матрицы. Устойчивость метода Гаусса. Опуская обременительные преобразования в методе обратного анализа ошибок округления, отметим, что возмущенная система метода Гаусса имеет вид A  F x  b . Запишем оценку нормы матрицы возмущения: F   1.01n 3  n 2 g A  A  . Вид этой оценки удовлетворял бы критерию устойчивости Уилкинсона, если бы множитель g(A) имел небольшое значение. Поясним смысл множителя g(A). Пусть A k  aijk обозначает матрицу, полученную из A после k шагов исключения. Обозначим   a0  max aij , ak  max aijk i,j . i,j Тогда g A   max 0 k  n 1 ak a0  . Следовательно, g(A) показывает, во сколько раз могут возрасти элементы матрицы A в ходе исключения переменных по сравнению с их исходным уровнем. По этой причине g(A) называют коэффициентом роста матрицы A . Элементы активной части матрицы Ak в методе Гаусса вычисляются по формуле a ijk  a ijk 1  a ikk 1 a kjk 1 a kkk 1 . Для ограничения роста элементов матрицы в процессе гауссова исключения желательно, чтобы поправочные члены ( k 1) ( k 1) ( k 1) aik a kj / a kk ( k 1) в этой формуле были не слишком большими. Это достигается процедурой выбора элемента akk , который называют главным. Выбор главного элемента по столбцу. В этом случае ограничение роста элементов матрицы Ak на k–м шаге гауссова исключения достигается перестановкой строк таким образом, чтобы гарантировать неравенство a ikk 1 a kkk 1  1 ; i  k,n . С этой целью при исключении переменной xk в качестве главного элемента выбирается элемент ( k 1) ark матрицы Ak-1 по правилу k 1 , a rkk 1  max a ik k i  n т. е. наибольший по модулю элемент в k–м столбце матрицы Ak-1 (рис. 3.1). Строки r и k переставляются и только после этого выполняется k–й шаг исключения прямого хода Гаусса. При столбцовой стратегии выбора главных элементов справедлива такая A k 1 оценка для значения параметра ak, определяющего коэффициент роста: k a k  2a k 1 . n  1 Она допускает, что an 1  2 a0 и, следовательно, коэффициент роста g A   2n1 . k Рис. 3.1. Матрица A в начале k-го шага 15 По этой причине метод Гаусса с выбором главного элемента по столбцам является условно устойчивым. Несмотря на это, он широко используется на практике, так как g(A) редко достигает своего верхнего предела. Выбор главного элемента по всей матрице. В этой стратегии в качестве главного элемента при исключении неизвестной xk выбирается элемент arl(k 1) по правилу ( k 1) a rl A k 1 k * r ( k 1)  max a ij , k i  n , k  j n т. е. наибольший по модулю элемент в квадратной k  k  подматрице матрицы Ak-1 (рис. 3.2). Строки k и r, а также столбцы k и l переставляются и далее выполняется k–й шаг исключения. Такая стратегия гарантирует выполнение неравенства ( k 1) aij(k 1) / akk  1; i  k,n; j  k,n и, следовательно, ограничивает рост элементов в процессе исключения Гаусса. Рис. 3.2 Область выбора главного Оценка коэффициента роста элементов матрицы A в этом случае элемента на k-м шаге исключения имеет более благоприятный вид: g A   n . Точность метода Гаусса. Привлекая оценку нормы матрицы возмущения, можно записать, что  x  1.01 n 3  n 2 g A  . Анализ неравенства позволяет определить пути повышения точности метода Гаусса: выбор главных элементов, работа с числами удвоенной длины, переобусловливание системы линейных алгебраических уравнений. k l Устойчивость и точность метода LU-факторизации. Метод LU-факторизации, если не делать никаких специальных усилий, характеризуется такими же оценками нормы матрицы возмущения и относительной ошибки, как и метод Гаусса. В то же время этот метод наряду с упоминавшимися при рассмотрении метода Гаусса направлениями повышения точности решения имеет еще одну возможность: использование так называемой операции скалярного накопления. Эта операция во многих ЭВМ реализована аппаратно либо моделируется специальными подпрограммами. Операция скалярного накопления привлекается при вычислении скалярного произведения n  k k . k 1 Пусть основным режимом работы ЭВМ является арифметика с плавающей точкой и tразрядной мантиссой машинного слова. Рассмотрим вычисление скалярного произведения. Произведение чисел βk и γk, выполненное арифметическим процессором, имеет 2t–разрядную мантиссу. Вместо того, чтобы округлить это произведение до t разрядов и прибавить результат к хранимой t–разрядной промежуточной сумме Sk-1, режим накопления предусматривает удержание 2t разрядов как для Sk-1, так и для βkk. Очередное суммирование S k  S k 1   k  k выполняется с 2tразрядными числами. Оно также сопровождается ошибкой округления, но эта ошибка затрагивает лишь последний разряд 2t–разрядного результата. Возвращение к t–разрядному представлению происходит лишь после вычисления скалярного произведения. Анализ расчетных соотношений метода LU-факторизации свидетельствует о том, что они реализуются с использованием операции скалярного накопления. Оценка нормы матрицы возмущений метода LU-факторизации при использовании операции скалярного накопления имеет вид F   g A  A  . Устойчивость вычислительной схемы метода определяется величиной коэффициента роста элементов матрицы A. Чтобы уменьшить величину g(A), в методе LU-факторизации используется 16 столбцовая процедура выбора главных элементов. Погрешность решения метода δ x  g (A) . Устойчивость и точность метода Холесского. Этот метод применяется для решения линейных систем с симметрической положительно определенной матрицей A. Можно показать, что в схеме Холесского коэффициент роста элементов матрицы A равен 1. Действительно, из соотношения ~2 ~2 ~2 a 2jj  l j1  l j 2    l jj ~ следует, что l ji  a jj . Если учесть, что в симметрической положительно определенной матрице A наибольший по модулю элемент положителен и находится на главной диагонали, нетрудно прийти к сформулированному заключению: g(A)=1. Таким образом, чтобы ошибки округления при реализации схемы Холесского были малыми, не требуется никаких перестановок строк либо столбцов. Наименьшая погрешность численного решения достигается при привлечении операции скалярного накопления. В этом случае F   A  и  x   . Прямые методы для систем с разреженными матрицами. Во многих приложениях приходится решать систему Ax  b с разреженной матрицей A. Назовем [n  n] –матрицу A разреженной, если число ее ненулевых элементов z  n 2 . При решении таких систем естественно хранить в памяти ЭВМ только ненулевые элементы. Можно надеяться, что в случае разреженных матриц как вычислительная работа, так и требуемый объем памяти значительно сократятся. Проблема состоит в том, что на этапе исключения неизвестных в методе Гаусса либо на этапе построения L и U матриц в компактной схеме Гаусса могут возникать Рис. 3.3. LU-факторизация диагональной матрицы со строчным и столбцовым окаймлением новые ненулевые элементы. Обратимся к следующему простому при-меру (рис. 3.3). Хотя исходная матрица A здесь имеет ничтожное число ненулевых членов, матрицы L и U имеют такую же структуру Рис. 3.4. LU-факторизация преобразованной диагональной матрицы со строчным и столбцовым окаймлением ненулевых элементов, как и для полностью заполненной матрицы. В 17 то же время перестановкой строк и столбцов систему можно преоб-разовать к виду, который сохраняет свойство разреженности исходной матрицы в матрицах разложения (рис. 3.4). В общем случае, конечно, не всегда удается полностью сохранить разреженность. Поэтому перестановки строк и столбцов выбирают такими, чтобы заполняемость матрицы в процессе исключения переменных была минимальной. На языке матричного описания эта задача выглядит следующим образом. Пусть P и Q – n  n матрицы перестановок строк и столбцов, причем QQ T  I . Матрицы P и Q 0 0 1 0 легко получить из единичной матрицы, в которой переставлены 0 1 0 0 Р= соответствующие строки (столбцы). Например, матрица перестановки 1 0 0 0 первой и третьей строк [4  4] -матрицы приведена на рис. 3.5. 0 0 0 1 Преобразуем систему линейных алге-браических уравнений Рис. 3.5. Пример матрицы перестановки следующим образом: PAQQ T x  Pb или By  d, где B  PAQ – переупорядоченная форма A, d  Pb , y  Q T x . К системе By  d будем применять прямой метод решения (его эффективность зависит от выбора матриц P и Q) и далее найдем x  Qy , так как Q T  Q 1 . Задача осложняется тем, что перестановка строк и столбцов определяющим образом влияет не только на разреженность матриц разложения, но и на численную устойчивость прямого метода. Симметричные матрицы. Если матрица A линейной системы является симметрической положительно определенной, то процесс исключения переменных по схеме Холесского является численно устойчивым при любом выборе главных элементов из числа диагональных членов. Таким образом, если положить Q  P T , т. е. при перестановке k–й и j–й строк осуществлять перестановку k–го и j–го столбцов, то выбор матрицы P можно производить только исходя из требования минимизации заполнения матрицы LT . Задача численного решения линейной системы сводится в этом случае к выполнению трех последовательных шагов: Шаг 1. Формирование матрицы B  PAP T и вектора d  Pb . Шаг 2. Решение упорядоченной системы By  d . Шаг 3. Вычисление вектора x  Qy . Рассмотрим более подробно, как реализуется Шаг 1. Вообще говоря, для [n  n] -матрицы A существует n ! различных упоря-дочиваний, из которых одно или более оказывается оптимальным в смысле заполнения. Задача отыскания оптимальных упорядочиваний является чрезвычайно трудной и практически нерешенной на сегодняшний день. Существующие процедуры реализации шага 1 являются эвристическими. Они обеспечивают понижение заполнения матрицы LT , однако не гарантируют достижение точного минимума. Опишем в общих чертах одну из таких процедур – алгоритм минимальной степени. Основная идея алгоритма минимальной степени заключается в локальной минимизации заполнения за счет выбора главного элемента в той строке и столбце, которые обеспечивают внесение наименьшего числа ненулевых элементов в множитель LT . Рассмотрим, например, k–й шаг исключения прямого хода Гаусса. В поддиагональных позициях столбцов с 1–го по (k-1)–й расположены нулевые элементы. Для исключения переменной xk на этом шаге k–я строка, нормированная соответствующим образом, вычитается из всех тех строк, которые имеют ненулевые элементы в k–м столбце. Следовательно, чтобы минимизировать число новых ненулевых элементов, достаточно просмотреть активную подматрицу матрицы Ak-1, выбрать строку с наименьшим числом ненулевых элементов, назовем ее l–й строкой, и переставить как строки, так и столбцы с номерами l, k. После такой перестановки в матрицу Ak-1 вносятся новые ненулевые элементы, т. е. формируется портрет матрицы Ak , и процесс имитации исключения переменных продолжается. Выполнив имитацию (n-1)–го шага исключения, построим в итоге из 18 единичной матрицы I искомую матрицу перестановок P. Матрицы общего вида. На k–м шаге алгоритма Гаусса при решении линейных систем с разреженной матрицей A произвольной структуры необходимо выбирать в качестве главного элемента такой элемент a (k 1) , который прежде всего обеспечивает устойчивость вычислительного l ,m процесса и в то же время по возможности минимизирует заполняемость матрицы Ak. Это достигается введением специального барьера – числа u≥1 и выделением множества элементов S k 1 из активной подматрицы матрицы Ak-1, которые удовлетворяют условию u al(,km1)  max ai(,kj1) , k  i  n, k  jn где al(,km1) – элемент множества Sk-1. Для каждого из таких элементов вводится числовая оценка M l , m, k  r (l,k )  1c(m,k )  1 , где r(l,k) – число ненулевых элементов в l–й строке активной части матрицы Ak-1, c(m,k) – число ненулевых элементов в m–м столбце этой же части матрицы Ak-1. В качестве главного элемента на k–м шаге прямого хода Гаусса выбирается такой элемент a *k,1*  S k 1 , l m для которого M l *, m*, k  min M l , m, k . k  l  n, k  mn Лекция 4 ИТЕРАЦИОННЫЕ МЕТОДЫ ЦЕЛЬ ЛЕКЦИИ: Дать общую схему линейных итерационных методов первого порядка и построить на её основе методы простой итерации, Гаусса–Зейделя, последовательной релаксации; сравнить эти методы по вычислительным затратам на итерации. В случае больших разреженных СЛАУ предпочтение отдаётся итерационным методам, т. к. они, во-первых, не приводят к появлению на итерациях новых ненулевых элементов, во-вторых, оказываются более эффективными по затратам машинного времени. Схема итерационных методов. Пусть требуется решить систему Ax  b ; b,x  R n ; A  R n n . Предположим, что Q – неособенная n  n -матрица, рассматриваемая как апроксимация матрицы A, для которой решение системы Qu  d не представляет особого труда. Это имеет место тогда, когда матрица Q является нижней треугольной, верхней треугольной, трёх-диагональной, произведением конечного числа таких простых матриц и т.п. Матрицу Q называют матрицей расщепления. Перепишем исходную систему в виде Qx  Q  A x  b . Построим итерационный процесс Qx k 1  Q  A x k  b ; k  0 ,1,2 , Разрешим его относительно x k 1 : x k 1  Q 1 (Q  A )x k  Q 1b , 19 или x k 1  (I  Q 1A )x k  Q 1b . Введем обозначения: B  I  Q1 A ; d  Q 1b . Тогда итерационную схему можно представить в виде x k 1  Bx k  d ; k  0 ,1,2 , Это линейная схема первого порядка. Поскольку итерационная матрица B на итерациях является постоянной, то такую схему называют стационарной. Достаточное условие сходимости этой итерационной схемы имеет вид B 1. Убедимся в этом. Обозначим вектор погрешности на k-м шаге через ε k , на (k+1)-ом шаге – ε k 1 , точное решение – x  . Тогда x k 1  x   ε k 1 ; x k  x   ε k . Учет этих соотношений в итерационной схеме позволяет записать x   ε k 1  B( x   ε k )  d , ε k 1  Bε k , ε k 1  B k 1ε 0 . Отсюда ε k 1  B При B 1 k 1 ε0 . видно, что ε k 1 k  0 , т. е. метод сходится к решению. Построим теперь  вычислительные схемы, выбирая конкретный способ формирования матрицы раcщепления. Метод простой итерации (Метод Якоби). Пусть требуется решить систему Ax  b ; x,b  R n . Представим A  A   D  A  , где D – диагональная матрица с диагональными членами матрицы A; A  – часть матрицы A, лежащая ниже центральной диагонали, A  – часть матрицы A, лежащая выше центральной диагонали. Тогда A  D  A x  b , или Dx  A   A  x  b . Запишем итерационный метод   Dx k 1   A   A  x k  b ; k  0,1,2, Разрешим его относительно x k 1 : x k 1  D 1 ( A   A  ) x k  D 1b ; k  0,1, Эта форма удобна для реализации метода Якоби. Здесь итерационная матрица B  D 1 A   A    D 1 A  D  I  D 1A  . Матрица расщепления Q  D . Следовательно, при построении метода Якоби в качестве матрицы ращепления используется диагональ матрицы СЛАУ. Запишем метод Якоби в координатной форме для системы общего вида: a11 x1  a12 x2    a1n xn  b1, a 21 x1  a 22 x2    a 2n xn  b2 ,  a n1 x1  a n 2 x2    ann xn  bn . Нетрудно убедиться, что метод Якоби в координатной форме x ik 1  i 1 1   bi   a ij x kj  j 1 a ii  n  j i 1   a ij x kj  ; i  1,n 20 есть не что иное, как разрешение каждого из уравнений системы относительно одной из компонент вектора. Из первого уравнения системы выражается x1 и его значение принимается за значение x1k 1 . Из второго уравнения определяется x2 и его значение принимается за x k2 1 и т. д. Переменные в правой части этих соотношений при этом полагаются равными их значениям на предыдущей итерации. Рассмотрим пример системы линейных алгебраических уравнений с разреженной (пятидиагональной) матрицей A (рис. 4.1) и оценим эффективность решения такой системы методом Якоби. Будем хранить матрицу A в виде пяти одномерных массивов g, f, a, b, c. Тогда решаемая система в терминах указанных массивов примет вид g i  m x i  m  f i 1 x i 1  a i x i  bi x i 1  c i x i  m  q i ; i  1, n . В этой записи следует учесть, что отдельные коэффициенты являются нулевыми. В частности g i  m  0, f i 1  0 при i=1; g i m  0 при i  2,m ; ci  0 при i  n  m  1,n  1 ; bi  0, ci  0 при i  n. m m i-я строка g i m f i 1 a i bi ci xi = x A qi q Рис. 4.1. Система с пятидиагональной матрицей Реализация метода Якоби для такой системы линейных алгебраических уравнений осуществляется просто: xik 1    1 qi  g i  m xik m  f i 1 xik1  bi xik1  ci xik m ; i  1,n; k  0,1,2, ai Видно, что для вычисления компоненты вектора решения необходимо выполнить четыре операции умножения и сложения и одну операцию деления. Полное время одной итерации метода Якоби оценивается соотношением T Я  (4t у 4t сл t дел )n. Метод Гаусса–Зейделя. Решаемую систему запишем в виде ( A   D  A  )x  b . Итерационная схема Гаусса–Зейделя также следует из этого пред-ставления системы: ( A   D)x k 1  b  A  x k ; k  0,1,2,... , или Dx k 1   A  x k 1  A  x k  b ; k  0,1,2,... . Последняя форма такого итерационного метода удобна для реа-лизации. Приведем метод Гаусса–Зейделя к стандартному виду: x k 1  ( A   D) 1 A  x k  ( A   D) 1 b ; k  0 ,1,... Стандартная форма метода позволяет выписать его итерационную матрицу и провести над ней очевидные преобразования: B  ( A   D) 1 A   ( A   D) 1 ( A  D  A  )  I  ( A   D) 1 A . Матрица расщепления Q  A  D cодержит всю нижнюю треугольную часть матрицы A. Запишем метод Гаусса–Зейделя в координатной форме для системы общего вида: 21 i 1 n   x ik 1   b i   a ij x kj 1   a ij x kj  a ii ; i  1,n; k  0,1,2,... j 1 j  i 1   Координатная форма метода Гаусса–Зейделя отличается от координатной формы метода Якоби лишь тем, что первая сумма в правой части итерационной формулы содержит компоненты вектора решения не на k-й, а на (k+1)-й итерации. Но именно это означает, что матрица расщепления метода Гаусса–Зейделя включает не только диагональные члены матрицы A, но и все члены, расположенные ниже главной диагонали. Оценим, как и ранее, временные затраты метода при решении тестовой задачи  СЛАУ с пятидиагональной матрицей. Реализация метода Гаусса–Зейделя для такой системы: xik 1  1 ai q  g i k 1 i  m xi  m   f i 1 xik11  bi xik1  ci xik m ; i  1,n; k  0,1,2,... позволяет констатировать, что временные затраты на итерации метода TГ  З  (4t у  4tсл  tдел )n оказываются такими же, как и в методе Якоби. Ключ к пониманию более высокой эффективности метода Гаусса–Зейделя  количество итераций, требующееся для получения решения с заданной точностью. Более высокая сходимость к решению в методе Гаусса–Зейделя достигается за счет выбора матрицы расщепления, лучше аппроксимирующей матрицу A. Метод последовательной релаксации. Воспользуемся методом Гаусса–Зейделя, записанным в форме Dx k 1   A  x k 1  A  x k  b , или x k 1 D 1 ( A  x k 1  A  x k  b) . Удобная для реализации итерационная схема релаксационного метода имеет вид x k 1  D 1 ( A  x k 1  A  x k  b)  1   x k ;k  0,1,2, , где ω – некоторый параметр. Преобразуем ее к стандартной форме:   x k 1  D  ω A   1 ωA  1  ωDx k   Dω  A   1 b. В свою очередь итерационная матрица  B  D  A    ωA  ωA 1  1 D  A   A . ω    D  I   В итоге матрица расщепления QD A  . При 1    2 указанная итерационная схема соответствует методу последовательной верхней релаксации, при 0    1 – методу последовательной нижней релаксации. Координатная форма метода последовательной релаксации: i 1 n   k x ik 1    bi   a ij x kj 1   a ij x kj  a ii  1  ωxi ; i  1,n; k  0 ,1,2 ,... j  1 j  1  1   Метод последовательной релаксации для системы с пятидиа-гональной матрицей: 1  k xik 1  ω qi  g i  m xikm1  f i 1 xik11  bi xik1  ci xik m   1  ωx i ,  ai    i  1,n; k  0 ,1,2 , Время выполнения итерации в этом случае Tп. р.  (6t у  5tсл  tдел )n . Метод касательных. Метод касательных (Ньютона) применяется для уточнения действительных корней нелинейного уравнения. 22 Имеем уравнение F(X) = 0. Если F(X)  D[A, B], и известно F(А) F(В) < 0 (F(А) и F(В) имеют разные знаки) и F11(X) не меняет знака на интервале [А,В], то уравнение F(X) = 0 имеет решение Х*  [A, B], F(X*)  0. Процесс уточнения действительного корня нелинейного уравнения методом касательных заключается в том, что строится касательная к кривой F(X) и определяется точка пересечения касательной с осью абсцисс, координата этой точки используется в качестве уточнения корня. Затем к кривой F(X) в точке последнего уточнения корня строят очередную касательную и определяется точка пересечения касательной с осью Х – очередное уточнение корня уравнения. Процесс уточнения корня (построение касательных) продолжается до тех пор пока два ближайших уточнения будут отличаться на величину не более , точность вычисления корня, | Х(i) – X(i-1) | <  где i - номер итерации.. Ограничение по применению метода. Метод касательных можно применять для уточнения действительного корня, если на интервале [А, В] функция F(X) удовлетворяет следующим свойствам: интервал [А,В] должен быть достаточно мал, чтобы на его длине график функции F(X) не имел горизонтальных участков, участков с малым наклоном касательных к кривой F(X) и экстремумов, функция F(X) должна быть монотонная, F11(А) F11(В) > 0; за начальное приближение Х0 принимается одна из границ интервала [А,В], где F(X) и F11(X) имеют одинаковые знаки, Х0 = А при условии F(А) F11(А) > 0 (Х0 = B при условии F(B) F11(B) > 0). Алгоритм метода. Исходные данные: начальное приближение корня Х, точность вычисления корня  . Организовать вычисление поправки к уточнению корня, Х = F(X) / F1(X), Уточнить значение корня. В качестве рекуррентной формулы метода касательных используют формулу вычисления координаты пересечения касательной с осью Х, Х = Х - Х, Проверить условие продолжения уточнения корня, | Х | >  . Если заданное условие принимает значение «истина», то необходимо продолжить уточнение корень уравнения с пункта 2. Если заданное условие принимает значение «ложно», то корень Х найден с заданной точностью , организовать вывод значения корня Х и прекратить уточнение корня. Часть 3. НЕЛИНЕЙНЫЕ АЛГЕБРАИЧЕСКИЕ УРАВНЕНИЯ И СИСТЕМЫ Лекция 5 ЧИСЛЕННОЕ РЕШЕНИЕ УРАВНЕНИЙ. МЕТОД ИТЕРАЦИЙ. ЦЕЛЬ ЛЕКЦИИ: Для уравнения f ( x)  0 ввести постановки задач отделения и уточнения корней, сформулировать лемму об оценке погрешности приближенного решения, построить метод простой итерации, получить для него оценку точности, доказать его сходимость. Пусть требуется решить уравнение f ( x )  0 , x  [a,b ] , т. е. найти все корни xi , удовлетворяющие этому уравнению на отрезке [a,b] . Задача численного решения уравнения f ( x)  0 сводится, во-первых, к отделению корней, вовторых, к последующему уточнению корней. Отделение корней. Отделить корни уравнения f ( x)  0, x  [a,b] значит заключить каждый корень xi в интервал 23 ai  xi  bi , для которого выполняются условия: 1. f (ai ) f (bi )  0 ; f ( bi ) 2. f (x) – знакопостоянная функция для x  [ai ,bi ] . Первое неравенство обеспечивает наличие в интервале [ai ,bi ] хоai тя бы одного корня, второе условие гарантирует единственность x bi x i корня (см. рис. 4.1). f (ai ) Для отделения корней можно использовать аналитический и табличный способы. Рис. 7.1. И ллю страция отд еления корня Аналитический предполагает исследование функции f ( x)  0, x  [a,b] методами математического анализа, последующее построение графика функции, из которого и определяются интервалы, содержащие единственный корень. Недостаток аналитического способа – неалгоритмизуемость. Табличный метод предполагает составление таблицы значений f ( xi ) , i  0,1,2, ,n , причем xi 1  xi  h, h  (b  a ) / n . Из таблицы на основании условий f ( x) f ( xi ) f ( xi 1 )  0, sign f ( xi )  sign f ( xi 1 ) алгоритмически определяются искомые интервалы [ai ,bi ] , ai  xi , bi  xi 1 . Чтобы не потерять корни, интервал отделения h должен быть достаточно малым. Уточнение корней. Пусть x  - корень уравнения f ( x)  0, x  [a,b] . В дальнейшем [a,b] – интервал, содержащий единственный корень. Уточнить корень с заданной точностью  значит найти приближенное значение xk  [a,b] такое, что x   xk   . Для решения сформулированной задачи необходимо уметь производить оценку абсолютной погрешности метода, т. е. величины x   xk . Лемма об оценке погрешности приближенного решения уравнения f ( x)  0 . Лемма. Пусть уравнение f ( x)  0 на отрезке [a,b] имеет корень x  . Пусть найдено некоторое его приближенное значение xk . Тогда x   xk  f ( xk ) m1 , где m1  min f ( x) . x[ a , b] Доказательство. Вычислим по теореме о конечных приращениях: f ( x  )  f ( xk )  f ( )( x   xk ) ,   [ x  ,xk ] . Очевидно, что f ( x  )  f ( xk )  f ( ) x   xk  m1 x   xk . Отсюда f ( xk )  m1 x   xk , или 24 x   xk  f ( xk ) m1 . Лемма доказана. f (x) f ( xk ) x xk xk x x   xk x x  xk Рис. 7.2. Влияние производной на погрешность приближенного решения Величину f ( xk ) называют невязкой. Из леммы следует, что судить о величине погрешности приближенного решения только по величине невязки нельзя. Необходимо учитывать и значение первой производной. Обратимся к рис. 4.2. Из этого рисунка следует, что одна и та же невязка приводит к существенно разным погрешностям приближенного решения, если производная в окрестности решения сильно отличается. Метод итераций. Для построения метода итераций преобразуем уравнение f ( x)  0, x  [a,b] к виду x   ( x) , x  [a,b] . Это можно сделать в общем случае так: x  x  hf (x) , или x   (x) , где  ( x)  x  hf ( x) . Постоянный множитель h выбирается при этом из условия сходимости метода (установим его позже). Пусть известно начальное приближение x0 . Тогда x1   ( x0 ) , x2   ( x1 ) ,  xk   ( xk 1 ). Приведенный способ построения числовой последовательности {xk } реализуется в методе итераций: xk   ( xk 1 ) , k  1,2,  . Рассмотрим, как ведет себя погрешность решения на итерациях метода. Обозначим xk 1  x    k 1, xk  x    k , где  k ,  k 1 - погрешности приближенного решения на двух соседних итерациях. Подставим представленные таким образом xk и xk 1 в итерационное правило: x    k 1   ( x    k ) . Разложим функцию  в ряд Тейлора в окрестности точки x  : x    k 1   ( x  )   k  ( x  )  O ( k2 ) . Пренебрегая остаточным членом O( k2 ) , получим соотношение, связывающее погрешность решения метода на двух соседних итерациях:  k 1   k  ( x  ) .  Сделаем некоторые выводы на основании этого приближенного равенства. Если  ( x )  1 , то можно ожидать, что  k 1   k и последовательность {xk } будет сходиться к решению, когда начальное приближение x0 выбрано достаточно близким к x  . 25  Если  ( x  )  1 , то скорее всего  k 1   k и метод будет расходиться, так как каждое после- дующее приближение будет отстоять от решения дальше, чем предыдущее.  При  ( x )  1 и  ( x  )  0 погрешности  k 1 и  k имеют одинаковые знаки. Сходимость будет монотонной.  При  ( x )  1 и  ( x  )  0 погрешности  k и  k 1 имеют разные знаки. Сходимость является немонотонной. Проиллюстрируем характер сходимости метода итераций графически. Образуем функции y1  x, y2   ( x) . В решении задачи значения этих функций совпадают. Первый пример (см. рис. 7.3,а) соответствует условиям  ( x  )  1,  ( x  )  0 и, следовательно, в этом случае наблюдается монотонная сходимость метода итераций. В следующем примере (см. рис. 4.3,б)  ( x  )  1,  ( x  )  0 , поэтому имеет место немонотонная сходимость. В последнем примере (см. рис. 4.3,в)  ( x  )  1 . При таких значениях производной метод итераций расходится. а) б) y1  x y1, y2 y1, y2 в) y2 (x) y2 (x)  x x2 x1 x0 x  x1 x3x y1, y2 y1  x x2 x0 x y1  x x2 y2 (x) x x1 x0 x Рис. 4.3. Иллюстрация сходимости метода итераций Теорема о сходимости и точности метода итераций. Теорема. Пусть уравнение f ( x)  0 приведено к виду x   (x) таким образом, что функция  (x) дифференцируема и выполняется условие  ( x)  q  1 для всех x  [a,b] . Тогда:   x ; последовательность {xk } k   ошибка x  xk  q xk  xk 1 1 q . Доказательство. Построим два соседних приближения в методе итераций: xk 1   ( xk ); xk   ( xk 1 ) . Очевидно, что xk 1  xk   ( xk )   ( xk 1 )   ( )( xk  xk 1 ) ,   [ xk ,xk 1 ] , или xk 1  xk   ( ) xk  xk 1  q xk  xk 1 . Таким образом, xk 1  xk  q xk  xk 1 . Запишем это неравенство для k=1,2,…,k в следующем виде: 26 x2  x1  q x1  x0 , x3  x2  q x2  x1  q 2 x1  x0 ,  (4.1) xk 1  xk  q k x1  x0 . Построим вспомогательный ряд x0  ( x1  x0 )  ( x2  x1 )    ( xk 1  xk )   (4.2) Частичная сумма членов этого ряда sk  xk , а значит {sk }  {xk } . По этой причине lim sk  lim xk . k  k  Последовательность частичных сумм ряда (4.2) совпадает с последовательностью приближений, вычисленных по методу итераций, а доказательство сходимости вспомогательного ряда эквивалентно доказательству сходимости метода итераций. Рассмотрим вспомогательный числовой ряд x0  x1  x0  q x1  x0  q 2 x1  x0    q k x1  x0   (4.3) Ряд (7.3) сходится абсолютно при q  1 . Но ряд (7.2) является мажорируемым к ряду (4.3) в силу неравенств (7.1). Следовательно, ряд (7.2) также сходится. Первая часть теоремы доказана. Напоминание. Ряд называется мажорируемым, если каждый его член по модулю не превосходит соответствующего члена некоторого сходящегося ряда с положительными членами. Получим теперь оценку погрешности приближенного решения. На основании леммы имеем: x   xk  f ( xk ) m1 . Найдем f ( xk ) . Так как f ( x)  x   ( x) , то f ( xk )  xk   ( xk ) . Теперь вычислим m1 . По определению m1  min f ( x)  min 1   ( x)  1  q . x[ a, b] x[ a ,b ] Тогда x   ( xk ) x  xk 1 q x   xk  k  k  xk  xk 1 1 q 1 q 1 q . Вторая часть теоремы также доказана. Замечание 1. Величину xk  xk 1 q /(1  q) можно использовать в качестве предельной абсолютной погрешности. Для ее расчета привлекаются два соседних приближения и минимальное по модулю значение первой производной. Предельная же абсолютная погрешность лежит в основе критерия завершения итерационного процесса. Критерий останова:  пред.   зад. , где  зад. – допустимая величина абсолютной ошибки, или xk  xk 1   зад. (1  q) / q . Замечание 2. Рассмотрим теперь, как выбирается h при построении функции  (x) . Напомним, что  ( x)  x  hf ( x) . Из теоремы следует, что для сходимости метода итераций необходимо, чтобы  ( x)  1, x  [a,b] . Значит 1  hf ( x)  1 или 1  1  hf ( x)  1 для x  [a,b] . Отсюда можно сделать вывод, что  при f x   0  при f x   0 0h 2 ; f x  2 h0. f x  Лекция 6 27 МЕТОД НЬЮТОНА ДЛЯ СИСТЕМ И ЕГО МОДИФИКАЦИИ ЦЕЛЬ ЛЕКЦИИ: Изучить метод Ньютона для решения систем нелинейных алгебраических уравнений; обсудить основные модификации метода: метод Ньютона с кусочно-постоянной матрицей, метод Ньютона–Рафсона, методы продолжения по параметру, метод Ньютона для плохо обусловленных систем, дискретный метод Ньютона. Предполагается знакомнство с методом Ньютона из курса информатики. Метод Ньютона для решения систем нелинейных уравнений. Пусть задана система нелинейных уравнений f1 ( x1,x2 , ,xn )  0, f 2 ( x1,x2 , ,xn )  0,  f n ( x1,x2 , ,xn )  0, решение которой достигается в точке ( x1 ,x2 , ,xn ) пространства R n . Обозначим F  ( f1,f 2 , ,f n ) T ; x  ( x1,x2 , ,xn ) T . Тогда исходная система запишется в виде F ( x)  0 , x  R n . Предположим, что известно k-е приближение x k к x  . Построим правило Ньютона вычисления (k+1)-го приближения в форме x k 1  x k  x k . Разложим функцию F(x k 1 )  F(x k  x k ) в ряд Тейлора в окрестности точки x k и сохраним в разложении два члена: F(x k 1 )  F(x k )  F (x k ) k x . x Полагая, что решение системы достигается на текущей итерации, относительно поправки x k получим систему линейных алгебраических уравнений: F(x k ) k k x  F ( x ) . x Тогда 1  F (x k )  k x     F(x ) ,  x  k и итерационное правило Ньютона решения системы нелинейных алгебраических уравнений запишется как 1 x k 1  F(x k )  k x   F(x ) ; k  0,1,2,  x   k Такой вид метода Ньютона неудобен на практике, потому что требует вычисления обратной матрицы, а эта операция достаточно трудоемка. На практике метод Ньютона реализуется в следующем виде: 1. Решается система линейных алгебраических уравнений и вычисляется вектор поправки: J ( x k ) x k   F ( x k ) , где J  F x – матрица Якоби системы; 2. Вычисляется (k+1)-е приближение x k 1  x k  x k , 3. Пункты 1, 2 повторяются для k=0,1,2,… до тех пор, пока не будет достигнута требуемая точность. Критерием завершения итерационного процесса служат условия 28 k k max f i   зад. , max xi   зад. , i1, n  i1, n  или в более общей форме F k   зад. , x k   зад. . Модификации метода Ньютона. Краткий обзор. Перечислим недостатки метода Ньютона, на устранение которых направлены различные его модификации:  трудность задания начального приближения, от которого метод сходится;  необходимость вычисления на каждой итерации матрицы Якоби, что может потребовать существенных вычислительных затрат;  необходимость решения на каждой итерации системы линейных алгебраических уравнений;  требование невырожденности матрицы Якоби. Рассмотрим модификации метода Ньютона, которые в той или иной мере устраняют перечисленные недостатки. Метод Ньютона с кусочно-постоянной матрицей. Чтобы уменьшить вычислительные затраты на итерации, матрица Якоби в этой модификации остается постоянной на протяжении нескольких шагов. Число шагов m, на которых J постоянна, задается в такой модификации в качестве параметра, либо момент перевычисления матрицы Якоби определяется условием F k 1   F k , в котором   0;1 , например   0.1 (матрица Якоби лишь при нарушении этого условия вычисляется заново). Эффективность метода достигается в этом случае не только путем сокращения числа расчетов матрицы Якоби, но главным образом за счет того, что на m итерациях метода приходится решать линейные системы с одной и той же матрицей. Метод Ньютона–Рафсона. Чтобы обеспечить сходимость метода от выбранного начального приближения x 0 , применяется модификация, называемая методом Ньютона–Рафсона. Вычисление (k+1)-го приближения в этой модификации осуществляется по правилу x k 1  x k   k x k , где  k – параметр, значение которого на k-й итерации выбирается из условия F k 1  F k . Стратегия выбора параметра  k на итерации может быть такой. Вначале принимается пробное значение  k  1 либо  k   k 1 и далее это значение видоизменяется до выполнения сформулированного условия. Это условие может потребовать многократного вычисления вектора F k 1 на текущей итерации. Очевидно, что при  k  1 метод Ньютона–Рафсона совпадает с методом Ньютона. Методы продолжения по параметру. Эти методы позволяют обеспечить сходимость метода Ньютона от выбранного начального приближения x 0 .Сущность методов продолжения по параметру заключается в замене исходной задачи последовательностью задач, каждая последующая задача при этом незначительно отличается от предыдущей. Последовательность строится таким образом, что первая система имеет решение x 0 , а последняя система совпадает с исходной задачей. Поскольку системы отличаются незначительно, то решение предыдущей задачи окажется хорошим начальным приближением для последующей. Решая такую последовательность задач методом Ньютона, получим в итоге решение исходной системы. Рассмотрим способ построения указанной последовательности задач. 29 Пусть при решении системы F ( x)  0 используется начальное приближение x . Заменим исходное уравнение F (x)  0 уравнением с параметром H (x ,t )  0 ; t  0,1 , которое при t  0 имеет решение x , а при t  1 совпадает с решением исходной задачи, т. е. H (x 0 ,0)  0; H (x ,1)  F (x  )  0 . В качестве H (x ,t ) можно выбрать функции H (x ,t )  (t  1)F(x 0 )  F(x) либо H (x ,t )  (x  x 0 )(1  t )  tF(x) . Разобьем отрезок 0,1 точками t0  0,t1,t 2 , ,t M  1 на M интервалов. Получим искомую последовательность систем: H (x ,ti )  0, i  0,1, ,M . Метод Ньютона для плохо обусловленных задач. Если матрица Якоби J  F x плохо обусловлена, погрешность решения линейной системы J ( x k ) x k   F ( x k ) может оказаться значительной из-за ошибок округления. Поэтому в случае плохо обусловленных матриц J при вычислении вектора поправок x k привлекают систему [ΘΕ  (1  Θ)J (x k )]x k  F(x k ) с числовым параметром Θ  0;1 , где E–единичная матрица. При Θ  0 модифицированная система совпадает с линейной системой стандартного метода Ньютона, при стремлении Θ к единице число обусловленности модифицированной системы также стремится к единице. Однако скорость сходимости соответствующего метода значительно ухудшается, поскольку при Θ  1 метод вырождается в метод простых итераций. Дискретный метод Ньютона. Дискретный метод Ньютона базируется на аппроксимации матрицы Якоби на основе вычисленных значений функции F(x) в ряде вспомогательных точек. Построим его. Как и прежде, будем использовать векторную запись решаемой системы уравнений F ( x)  0 . Пусть известно k-е приближение x k к решению x  . Аппроксимируем функцию F(x) линейной функцией: l (x)  Bx  b . Для численного определения матрицы B и вектора b потребуем, чтобы значения функций F(x) и l (x) совпадали в (n+1) вспомогательных точках x k,j , j  0,1,2, ,n , т. е. чтобы выполнялось равенство Bx k,j  b  F(x k,j ); j  0,n . Вычитая из первого равенства все последующие, получим соотношения B(x k,0  x k,j )  F (x k,0 )  F(x k,j ); j  1,n или в матричной форме BX  Ф , где [n  n] матрицы X и Ф имеют вид X  [x k,0  x k,1, x k,0  x k,2 ,  ,x k,0  x k,n ]; Φ  [F(x k,0 )  F (x k,1 ) , F(x k,0 )  F(x k,2 ) ,  ,F(x k,0 )  F(x k,n )]. 30 Следовательно, B  ΦX 1 . Условие l(x k,0 )  F(x k,0 ) позволяет найти вектор b : b  F(x k,0 )  ΦX 1x k,0 . Перепишем функцию l(x) , подставив найденные соотношения для B и b : l (x)  ΦX 1 (x  x k ,0 )  F(x k ,0 ). Примем x k  x k,0 . Будем искать x k 1 из уравнения l (x k 1 )  ΦX 1 (x k 1  x k )  F(x k )  0. Получим итерационную формулу дискретного метода Ньютона: x k 1  x k  XΦ 1F(x k ); k  0 ,1,2, . Заменяя обращение матрицы решением линейной системы, придем к реализуемому на практике алгоритму дискретного метода Ньютона: 1. Вычисляется вектор F k  F(x k ) , матрицы X k  X(x k ) и Φ k  Φ(x k ) . 2. Решается система линейных алгебраических уравнений Φ k v k  F k . 3. Вычисляется вектор поправки w k  Xk v k . 4. Вычисляется (k+1)-е приближение x k 1  x k  w k . 5. Пункты 1÷4 повторяются для k=0,1,2,… до получения решения с требуемой точностью. Применение дискретного метода Ньютона предполагает хранение [n  n] -матриц X и Ф . Однако на практике в качестве вектора x k,j выбирается вектор x k,j  x k,0  he j , где h – диагональная матрица параметров дискретизации, e j – j-й столбец единичной матрицы. Элементы h j матрицы h вычисляют по правилу h j  Mx kj ,0 , где M – константа (например, 0.1). 31 Часть 4. ИНТЕРПОЛИРОВАНИЕ ФУНКЦИЙ Лекция 7 ПОЛИНОМИАЛЬНАЯ ИНТЕРПОЛЯЦИЯ ЦЕЛЬ ЛЕКЦИИ: Пояснить сущность задачи интерполяции; дать общую схему ее решения; построить интерполяционный полином Лагранжа, ввести разделенные разности и получить на их основе интерполяционный многочлен Ньютона. Постановка задачи. Важным элементом численного анализа является задача интерполирования функций, которую приходится решать в следующих приложениях: машинная графика, автоматизация проектирования, обработка экспериментальных данных, управление, передача данных и др. При этом часто возникает задача восстановления функции на отрезке, если известны ее значения в конечном числе точек отрезка. Эти значения определяются на основании экспериментальных измерений либо в результате вычислений. Сформулируем математическую постановку задачи интерполяции. Пусть на отрезке [a,b] в точках a  x0  x1  x2    xn  b известны значения функции y(x) , равные y ( x0 )  y0 , y ( x1 )  y1,  , y ( xn )  yn . Требуется построить интерполянту – интерполяционную функцию f ( x) , x  [a,b] , совпадающую с функцией y (x) в точках x0 ,x1, ,xn : f ( xi )  yi ; i  0,n (7.1) и в то же время аппроксимирующую ее, естественно приближенно, на всем отрезке [a,b] . Основной вопрос интерполяции: как выбрать интерполянту f (x) и как оценить погрешность y ( x)  f ( x) ? Интерполяционные функции строятся, как правило, в виде линейной комбинации некоторых линейно-независимых элементарных функций Φk ( x) , k  0,n : n f ( x)   ck Φk ( x) , (7.2) k 0 где ck – коэффициенты, которые можно определить, используя условие (7.1). Из этого условия n  ck Φk ( xi )  yi ; i  0,n . (7.3) k 0 Относительно коэффициентов ck получили систему линейных алгебраических уравнений, порядок которой равен n+1. В силу линейной независимости элементарных функций Φk (x) определитель этой системы отличен от нуля, и решение для ck , k  0,n единственно. Таким образом, вычисляя ck из (7.3) и подставляя их в (7.2), можем получить значение f (x) в любой точке x  [a,b] . В качестве системы линейно-независимых функций Φk (x) чаще всего выбирают степенные функции, например Φk ( x)  x k . В этом случае f ( x)  Pn ( x) – полином степени n. Интерполяционная формула Лагранжа. При построении интерполяционной формулы Лагранжа в качестве Φk (x) используются полиномы Лагранжа lk (x) степени n, удовлетворяющие условиям 1, если k  i, lk ( xi )   i,k  0,n . 0, если k  i, Для выполнения второго условия полином степени n должен иметь вид 32 lk ( x)   ( x  x0 )( x  x1 )  ( x  xk 1 )( x  xk 1 )  ( x  xn ) , т. е. его корнями являются все узлы интерполяции, кроме k-го. Коэффициент  определим, используя первое условие: lk ( xk )   ( xk  x0 )( xk  x1 )  ( xk  xk 1 )( xk  xk 1 )  ( xk  xn )  1 . Отсюда находим  1 ( xk  x0 )  ( xk  xk 1 )( xk  xk 1 )  ( xk  xn ) и lk ( x )  n xx ( x  x0 )  ( x  xk 1 )( x  xk 1 )  ( x  xn ) i .  ( xk  x0 )  ( xk  xk 1 )( xk  xk 1 )  ( xk  xn ) i  0 xk  xi ik Решение системы (7.3) при использовании в качестве элементарных функций полиномов Лагранжа lk (x) имеет вид ck  yk . Таким образом, интерполяционный многочлен Лагранжа представится как n Pn ( x)  n n x  xi . i  0 xk  xi  yk lk ( x)   yk  k 0 k 0 ik Преобразуем его к виду, используемому на практике при вычислении значений функции:      n  n  y 1 k . Pn ( x)    ( x  xi )   n  i  0  k  0   ( x  x ) x  xk  k i    ii  0k    Интерполяционный многочлен совпадает с интерполируемой функцией только в точках В остальных точках имеет место погрешность интерполяции R( x)  y ( x)  Pn ( x) , которая оценивается величиной R( x)  y ( n 1) ( ) H ( x) /(n  1)!;   [a,b] , x0 ,x1,x2 , ,xn . n где H ( x)   ( x  xi ) . (Вывод оценки погрешности опустим). Погрешность интерполяции зависит от i 0 числа узлов интерполяции n, а также от их расположения на отрезке [a,b] . Наилучшими узлами интерполяции следует признать те xi , для которых max H ( x) принимает наименьшее значение. Недостатком интерполяционной формулы Лагранжа является то, что каждое слагаемое зависит от всех узлов интерполяции. При добавлении узла интерполяции и, следовательно, повышении порядка полинома необходимо вычислять не только слагаемое, относящееся к новому узлу, но и перевычислять заново слагаемые, относящиеся ко всем узлам интерполяции. Понятие конечных и разделенных разностей. Пусть на равномерной сетке в узлах x0 ,x1,x2 , ,xn заданы значения функции y(x) , т. е. yi  y ( xi ) , i  0,n; xi  xi 1  h; i  1,n . Конечными разностями первого порядка называют величины yi  yi 1  yi , i  0 ,n  1 . Разности конечных разностей первого порядка называют конечными разностями второго порядка: 2  yi  yi 1  yi , i  0 ,n  2 . Аналогично определяют конечные разности следующих порядков: k k 1 yi 1  k 1 yi , i  0,n  k .  yi   Схему вычисления конечных разностей удобно представить в виде следующего треугольника: 33 x0 x1 x2 x3 x4 y0 y1 y2 y3 y4 y0 y1 y 2 y 3 2 y0  2 y1 2 y2 3 y 0  3 y1 4 y0 Если узлы интерполяции – произвольные точки, то вводится понятие разделенных разностей. Разделенная разность первого порядка: y ( xi ,x j )  y ( xi )  y ( x j ) xi  x j . Разделенная разность второго порядка: y ( xi ,x j ,xk )  y ( xi ,x j )  y ( x j ,xk ) xi  xk . Разделенная разность k-го порядка: y ( xi ,xi 1, ,xi  k )  y ( xi , ,xr 1,xr 1, ,xi  k )  y ( xi , ,xs 1,xs 1, ,xi  k ) . Отметим xs  xr два важных свойства разде- ленных разностей: 1. Разделенная разность является симметричной функцией своих аргументов, т. е. значение разделенной разности не зависит от порядка перечисления узлов интерполяции, а зависит только от их значений; 2. Разделенная разность вычисляется и тогда, когда две или более точки интерполяции одинаковы. Разделенная разность первого порядка в случае, когда одинаковы две точки, имеет вид y ( xi ,xi )  lim y ( xi )  y ( x j ) xi  x j xi  x j  y( xi ) . Если точка xi среди узлов интерполяции встречается (m+1) раз, то разделенная разность m-го порядка на этих узлах y ( xi , ,xi )  y ( m) ( xi ) m! . Приведем теперь схему вычисления разделенных разностей x0 x1 x2 x3 x4 y0 y ( x0 ,x1 ) y1 y ( x1,x2 ) y2 y ( x 2 ,x3 ) y3 y ( x3 ,x4 ) y4 y ( x0 ,x1,x 2 ) y ( x1,x2 ,x3 ) y ( x2 ,x3 ,x4 ) y ( x0 ,x1,x2 ,x3 ) y ( x1,x 2 ,x3 ,x 4 ) y ( x0 ,x1,x 2 ,x3 ,x 4 ) Интерполяционная формула Ньютона. Пусть x  [a,b] . Построим первую разделенную разность y ( x,x0 )  y ( x)  y ( x0 ) x  x0 . Откуда находим y ( x)  y ( x0 )  ( x  x0 ) y ( x,x0 ) . Построим вторую разделенную разность y ( x,x0 ,x1 )  y ( x,x0 )  y ( x0 ,x1 ) x  x1 и выразим из нее y( x,x0 ) : y ( x,x0 )  y ( x0 ,x1 )  ( x  x1 ) y ( x,x0 ,x1 ) . Подставим y ( x,x0 ) в выражение для y (x) : y ( x)  y ( x0 )  ( x  x0 ) y ( x0 ,x1 )  ( x  x0 )( x  x1 ) y ( x,x0 ,x1 ) . Аналогично привлекаем следующие узлы интерполяции для построения интерполяционной функ34 ции. После использования всех узлов интерполяции: y ( x)  y ( x0 )  ( x  x0 ) y ( x0 ,x1 )  ( x  x0 )( x  x1 ) y ( x0 ,x1,x2 )     ( x  x0 )  ( x  xn 1 ) y ( x0 , ,xn )  ( x  x0 )  ( x  xn ) y ( x,x0 , ,xn )  Таким образом, интерполяционный мно-  Pn ( x)  Rn ( x). гочлен Ньютона имеет вид Pn ( x)  y0  n 1 i  i  0 k 0    y ( x0 , ,xi 1 )  ( x  xk ) . Непосредственной проверкой нетрудно убедиться, что yk  Pn ( xk ) , k  0,n . Погрешность интерполяции n Rn ( x)  y ( x,x0 , ,xn )  ( x  xi ) . i 0 Более эффективное вычисление значения функции по интерполяционной формуле Ньютона можно получить, если преобразовать ее к такому виду: y ( x )  y ( x0 )  ( x  x0 )[ y ( x0 ,x1 )  ( x  x1 )[ y ( x0 ,x1,x2 )     ( x  xn  2 )[ y ( x0 , ,xn 1 )  ( x  xn 1 ) y ( x0 , ,xn )]]. Интерполяционная формула Ньютона позволяет легко наращивать число узлов интерполяции, требуя при этом вычисления лишь дополнительных слагаемых. Например, добавление узла xn 1 приведет к вычислению слагаемого ( x  x0 )( x  x1 )  ( x  xn ) y ( x0 ,x1, ,xn ) . Лекция 8 СПЛАЙН-ИНТЕРПОЛЯЦИЯ ЦЕЛЬ ЛЕКЦИИ: Показать ограничения полиномиальной интерполяции, ввести понятие mсплайна, построить вычислительную схему интерполяции на основе кубического сплайна. Ограничения многочленной интерполяции. Многочленная интерполяционная функция очень чувствительна к выбору узлов интерполяции. Рассмотрим интерполяционную формулу Лагранжа n Pn ( x)   yk lk ( x); x  [a,b] , k 0 где n lk ( x)  x  xi . i  0 xk  xi  ik Оценим норму функции Pn (x) , которую определим, как Pn ( x)  max Pn ( x) . a  xb С этой целью выполним очевидные преобразования: Pn ( x)  n n k 0 k 0 n  yk lk ( x)   yk lk ( x)  (max yk )  lk ( x) . k k 0 Введем функцию Лебега: n n ( x )   l k ( x ) . k 0 35 Поскольку max y k  y ( x) , получаем оценку: k Pn ( x)  y ( x) n ( x) . Величина нормы функции n (x) зависит от распределения узлов интерполяции. Рассмотрим два случая. Случай 1. Узлы интерполяции xi на отрезке [a,b] распределены равномерно. В этом случае n ( x)  const  e n 2 (доказательство опускаем). При увеличении n многочлен Pn (x) может полностью оказаться непригодным для аппроксимации y (x) , т. к. Pn (x) неограниченно возрастает. Рассмотрим такой пример. Будем интерполировать функцию 1 y ( x)  1  25 x 2 , x  [1,1] полиномом n-го порядка, выбирая узлы интерполяции равномерно распределенными. Введем норму ошибки интерполяции: en  max y ( x)  Pn ( x) 1 x 1 и исследуем ее зависимость от порядка интерполирующего полинома. Для этого обратимся к численным результатам (см. табл.8.1). Таблица 8.1 n en n en n en 2 4 6 0.96 0.71 0.43 8 10 12 0.25 0.30 0.56 14 16 18 1.07 2.10 4.21 Видно, что с увеличением n погрешность интерполяции уменьшается вплоть до n  8 . Дальнейшее увеличение порядка полинома приводит к возрастанию en . Случай 2. Узлы интерполяции являются нулями полинома Чебышева. Их нетрудно получить, если на отрезке [a,b] построить полуокружность, разделить ее на n  1 равные части и спроектировать на отрезок середину каждой из них (рис. 8.1). В этом случае n ( x)  2 ln n  4 (доказательство опускаем). Обратимся к рассмотренному выше примеру (см. табл.8.2) Погрешность интерполяции теперь при возрастании порядка полинома монотонно падает. a b Рис. 8.1. Расположение нулей полинома Чебышева Таблица 8.2 n en n en n en 2 4 6 0.93 0.75 0.56 8 10 12 0.39 0.27 0.18 14 16 18 0.12 0.08 0.06 Многочленная интерполяция в точках Чебышева позволяет увеличивать точность приближения функции посредством увеличения порядка полинома. Однако привлекать узлы Чебышева не всегда удается. (Например, в узле Чебышева функция имеет особенность.) Чтобы избежать зависимости точности аппроксимации от локальных свойств функции, переходят к кусочнополиномиальной аппроксимации. 36 Сплайн-интерполяция. Назовем m-сплайнами полиномы невысокого порядка m, которыми аппроксимируется функция y(x) на интервалах [ xi 1,xi ] , i  1,n и которые сшиваются в узлах интерполяции на основе требования непрерывности интерполирующей функции и ее первых (m  1) производных. Рассмотрим наиболее широко используемую на практике кубическую сплайн-интерполяцию. В ней искомая функция y (x) на интервале xi 1  x  xi аппроксимируется полиномом третьей степени f ( x)  ai  bi ( x  xi 1 )  ci ( x  xi 1 ) 2  d i ( x  xi 1 )3 , x  [ xi 1,xi ]. (8.1) Аналогичные полиномы можно записать для всех n интервалов, т. е. соотношение (8.1) справедливо для i  1,n . Условия сшивки сплайнов в узлах интерполяции: f ( xi )  yi ; i  0,n , (8.2) f ( xi  0)  f ( xi  0); i  1,n  1 , (8.3) f ( xi  0)  f ( xi  0); i  1,n  1 . (8.4) Кроме того, необходимо задать граничные условия в точках x0 и xn . Это можно сделать несколькими способами. Здесь будем считать, что f ( x0 )  0; f ( xn )  0 , т. е. будем рассматривать интерполяцию так называемыми естественными сплайнами. Условие (8.2) приводит к уравнениям: f ( xi 1 )  ai  yi 1, i  1,n ; f ( xi )  ai  bi hi  ci hi2  d i hi3  yi , i  1,n ; hi  xi  xi 1, i  1,n . (8.5) Получили 2n уравнений относительно 4n неизвестных ai ,bi ,ci ,di , i  1,n . Вычислим производные: f ( x)  bi  2ci ( x  xi 1 )  3d i ( x  xi 1 ) 2 ; f ( x)  2ci  6di ( x  xi 1 ) . Построим уравнения, к которым приводит условие непрерывности первой производной в узлах i  1,n  1 . При кубической сплайн-интерполяции выражения для первой производной на соседних интервалах интерполирования имеют вид f ( x)  bi  2ci ( x  xi 1 )  3di ( x  xi 1 ) 2 , x  [ xi 1,xi ]; f ( x)  bi 1  2ci 1 ( x  xi )  3d i 1 ( x  xi ) 2 , x  [ xi ,xi 1 ]. Требование непрерывности первой производной в узлах i  1,n  1 (условие 8.3) приводит к уравнениям bi  2ci hi  3d i hi2  bi 1; i  1,n  1 . (8.6) В свою очередь, непрерывность второй производной в этих же узлах (условие (8.4)) позволяет записать 2ci  6d i hi  2ci 1; i  1,n  1 , или ci  3di hi  ci 1; i  1,n  1 . (8.7)     Наконец, из граничных условий f (0)  0 , f ( xn )  0 получим, что c1  0; cn  3d n hn  0 . (8.8) Соотношения (8.5), (8.6), (8.7), (8.8) составляют систему линейных алгебраических уравнений, всего 4n уравнений, относительно 4n неизвестных. Построим эффективный метод решения этой системы. Выразим из (8.7) d i : di  ci 1  ci ; i  1,n  1 , 3hi (8.9) 37 и подставим в (8.5): ai  hi bi  hi2 ci  hi3 (ci 1  ci )  yi ; i  1,n  1 . 3hi Учтем, что ai  yi 1, i  1,n  1 : yi 1  hi bi  hi2 ci  hi2 (ci 1  ci )  yi ; i  1,n  1 . 3 Выразим из этого соотношения bi : bi  yi  yi 1 hi  (ci 1  2ci ); i  1,n  1 . hi 3 (8.10) Подставим теперь bi ,bi 1,d i в формулу (8.6): yi  yi 1 hi c c  (ci 1  2ci )  2hi ci  3hi2 i 1 i  hi 3 3hi  yi 1  yi hi 1  (ci  2  2ci 1 ); i  1,n  1 . hi 1 3 Выполним очевидные преобразования и запишем результирующую систему в виде: y y y  yi 1  hi ci  2(hi  hi 1 )ci 1  hi 1ci  2  3 i 1 i  i ; h hi  i 1  (8.11) i  1,n  1 . Получили систему (n  1) линейных алгебраических уравнений относительно (n  1) неизвестных c1,c2 , ,cn 1 . Однако, если учесть, что из граничных условий c1  0, cn 1  cn  3d n hn  0 , то приходим к системе линейных алгебраических уравнений с треугольной симметрической матрицей a12c2  a13c3  b1, a22 c2  a23c3  a24c4  b2 , a33c3  a34c4  a35c5  b3,  an  2, n  2 cn  2  an  2, n 1cn 1  an  2, n cn  bn  2 , an 1, n 1cn 1  an 1, n cn  bn 1. Обратите внимание: матрица системы обладает диагональным преобладанием! Эту систему можно эффективно решать методом LLT -факторизации. Вычислив коэффициенты c2 , ,cn , из соотношений (8.9) находим di  ci 1  ci ; i  1,n  1 . 3hi Из условия (8.8) определяем d n : d n  cn / 3hn . Коэффициенты bi рассчитаем из (8.10): bi  yi  yi 1 hi  (ci 1  2ci ); i  1,n  1, hi 3 bn  yn  yn 1 2  hn cn , hn 3 так как cn 1  0 . Коэффициенты ai , i  1,n известны из (8.5). В результате для всех n кубических сплайнов определены коэффициенты ai ,bi ,ci ,d i . Расчет значений функции f (x) методом сплайн-интерполяции осуществляется следующим образом: 1. По описанной выше методике вычисляются коэффициенты ai ,bi ,ci ,d i ; i  1,n кубических сплайнов; 38 2. Находится интервал [ xi 1,xi ] , которому принадлежит данное x . Значение функции f (x) на этом интервале вычисляется из кубического сплайна f ( x)  ai  bi ( x  xi 1 )  ci ( x  xi 1 ) 2  di ( x  xi 1 )3 с параметрами ai ,bi ,ci ,di для интервала [ xi 1,xi ] . 39 Часть 5. ОПРЕДЕЛЕННЫЕ ИНТЕГРАЛЫ Лекция 9 ЧИСЛЕННОЕ РЕШЕНИЕ СИСТЕМ ОБЫКНОВЕННЫХ ДИФФЕРЕНЦИАЛЬНЫХ УРАВНЕНИЙ ЦЕЛЬ ЛЕКЦИИ: Дать постановку задачи численного решения системы обыкновенных дифференциальных уравнений с начальными условиями; провести классификацию численных методов по виду разностной схемы; ввести понятия локальной и интегральной точности, устойчивости по отношению к шагу интегрирования; построить простейший явный одношаговый метод, оценить его локальную погрешность, устойчивость к шагу интегрирования; пояснить свойство жесткости дифференциальных уравнений и непригодность для их решения классических явных методов. Общая характеристика численных методов. Пусть требуется решить систему обыкновенных дифференциальных уравнений (СОДУ) du1  f1 (u1,u 2 , ,u n ,t ) , dt  du n  f n (u1,u 2 , ,u n ,t ) dt с начальными условиями u1 (t0 )  u10 , u2 (t0 )  u20 ,  , un (t0 )  un0 . Это значит, что необходимо найти функции u1 (t ) , u2 (t ) ,  , un (t ) на отрезке t0  t  T такие, что u1 (t0 )  u10 , ,u n (t0 )  u n0 . Запишем сформулированную задачу, ее называют задачей Коши для системы обыкновенных дифференциальных уравнений, в векторном виде: du  f (u ,t ) , dt u(t0 )  u 0 . Решение задачи Коши заключается в нахождении такой траектории u(t ) , которая бы удовлетворяла условию u(t0 )  u 0 . Применение численных методов для решения задачи Коши предполагает приближенное вычисление значений u k  u(tk ) в точках t0 ,t1,t 2 , ,t n  T отрезка [t0 ,T ] . С этой целью СОДУ заменяют разностной схемой, из которой рекуррентно вычисляют приближенные значения u1,u 2 , ,u n . Запишем общий вид разностной схемы: y k 1   k (y k 1,y k , ,y k  p 1 )  0 . Она связывает искомое решение y k 1 в текущий момент времени tk 1 с построенными к этому моменту времени решениями y k , ,y k  p 1 в предыдущие моменты времени. Здесь y j – приближенное значение u j , т. е. y 0  u 0 , y1  u1, y 2  u 2 ,  . Выбор функции  k в разностной схеме определяет соответствующий численный метод. Приведем общую характеристику методов.  Если p  1 , т. е. y k 1   k (y k 1,y k )  0 , то метод является одношаговым. Он связывает решение в последующий момент времени с решением в предыдущий момент времени. Если p  1 , то метод является многошаговым.  Метод является явным, если функция  k не зависит от y k 1 и неявным в противном случае. В случае явного одношагового метода искомое решение y k 1 на текущем шаге интегрирования вычисляется тривиально: 40 y k 1   k (y k ) . Неявный одношаговый метод требует для расчета y k 1 решать систему в общем случае нелинейных алгебраических уравнений  (y k 1 )  y k 1   k (y k 1,y k )  0 , привлекая, например, итерационную процедуру Ньютона.  Различают локальную R и интегральную  погрешности метода. Локальная погрешность – ошибка на шаге интегрирования, интегральная погрешность – полная погрешность в текущий момент времени. Обратимся к графической иллюстрации погрешности. Пусть численно решается уравнение du  f (u,t ) , u (t )  u . dt В результате выполнения первого шага интегрирования (см. рис. 10.1) приближенное значение искомой функции в момент времени t1 равно u u0 кривая 1 R1=δ1 кривая 2 кривая 3 R2 t0 t1 t2 δ2 t Рис. 10.1. Иллюстрация локальной и интегральной погрешностей y1 . Это значение принадлежит некоторой интегральной кривой 2, являющейся решением дифференциального уравнения при другом, отличном от u0 , начальном условии. Разница между приближенным значением y1 и точным u (t1 ) составляет локальную погрешность. Очевидно, что на первом шаге локальная погрешность совпадает с интегральной погрешностью. Второй шаг интегрирования приводит к приближенному значению y2 , которое принадлежит интегральной кривой 3. Видно, что локальная погрешность на этом шаге, равная разности y2 и значения ординаты кривой 2 в момент времени t2 , существенно отличается от интегральной погрешности – разности между приближенным y2 и точным u (t2 ) значениями. Наиболее объективной характеристикой точности метода является величина интегральной погрешности. К сожалению, выполнить ее оценку крайне сложно. На практике при выборе шага интегрирования обычно используют локальную погрешность численного метода.  Если величина шага интегрирования  k  t k 1  t k ограничивается только допустимой локальной погрешностью R( k ) метода, то такой метод является абсолютно устойчивым к шагу интегрирования. Метод, обладающий таким свойством, позволяет осуществлять выбор шага интегрирования, исходя лишь из требуемой точности. Условно устойчивый метод характеризуется тем, что шаг интегрирования ограничивается не только допустимой локальной погрешностью, но и определенными свойствами решаемой системы. Последние ограничения являются весьма обременительными для некоторых классов задач. Рассмотрим разностные схемы наиболее широко используемых на практике методов численного интегрирования. Явный метод Эйлера. Рассмотрим u k 1  u(tk 1 ) , где tk 1  tk   k ,  k – текущий шаг интегрирования. Разложим функцию u k 1  u(tk   k ) в ряд Тейлора в окрестности точки tk : 41 u k 1  u k   k  2 d 2u  k3 d 3u du  k   . dt t  t k 2 ! dt 2 t  t k 3! dt 3 t t k Ограничившись в этом разложении двумя членами, получим разност-ную схему метода Эйлера y k 1  y k   k f (y k ,t k ) , k  0,1,2, . Локальная погрешность метода Эйлера составляет величину R ( k )   k2 d 2 u 2 dt 2 t t k  O( k3 ) . В вычислительной математике численные методы решения обык-новенных дифференциальных уравнений принято характеризовать порядком точности. Определение. Если локальная погрешность численного метода R ( k )  O( kp 1 ) ,то порядок точности такого метода равен p . Метод Эйлера является методом первого порядка. Приведем геометрическую интерпретацию явного метода Эйлера для задачи Коши du  f (u,t ) , u (t )  u dt (см. рис. 10.2). Приращение на шаге интегрирования – катет прямоугольного треугольника, лежащий против угла, тангенс которого равен значению производной в предыдущий момент времени. Вторым катетом этого треугольника является текущий шаг интегрирования. u y1 R( ) f ( y 0 , t0 ) u0 y0 t0 τ t1 t Рис. 10.2. Геометрическая иллюстрация явного метода Эйлера Оценим устойчивость метода Эйлера по отношению к шагу интегрирования. Для этого рассмотрим линейную автономную систему du  Au dt с отрицательно определенной [n  n] матрицей A простой структуры. Отрицательная определен- ность матрицы означает, что все собствен-ные значения i матрицы действительны и отрицатель 0 . ны, т. е. i  0 , i  1,n . В этом случае все решения u(t,u 0 ) t  Применим для решения этой системы метод Эйлера с постоянным шагом  : y k 1  y k  Ay k  (Ε  A )y k , k  0 ,1,2, . Здесь E – единичная матрица соответствующей размерности. Из алгебры известно, что для любой неособенной матрицы A простой структуры существует такая неособенная матрица S , которая преобразованием подобия приводит матрицу A к диагональному виду: S 1AS  diag {i } . Преобразуем вычислительную схему метода Эйлера следующим образом: S 1y k 1  S 1 (Ε  A)SS 1y k . Введем замену переменных z  S 1y . Тогда z k 1  S 1 (Ε  A)Sz k , или z k 1  [Ε   diag{λi }]z k , k  0 ,1,2 , . 42 Запишем это соотношение для i-й компоненты вектора z : z ki 1  (1  i ) z ki , k  0,1,2, , или z ki 1  (1  i ) k 1 z0i , k  0,1,2, , где z 0  S 1y 0 , т. е. определяется начальным условием. Нетрудно видеть, что z ki 1  0 , k  если 1  i  1 . Именно этим свойством обладает решение автономной системы с отрицательно определенной матрицей. Отсюда приходим к требованиям 1  1  i  1 , при этом неравенство 1  i  1 приводит к естественному условию   0 , т. к. i  0 , а неравенство 1  i  1  к условию   2 / i .  0 , необходимо при выборе шага интегрирования выполнить условие Очевидно, чтобы z k 1 k   2 . max{i } i Таким образом, явный метод Эйлера по отношению к шагу интегрирования является условно устойчивым. Явление жесткости. Ограниченная устойчивость численного метода является серьезным недостатком при решении так называемых жестких систем. Рассмотрим дифференциальное уравнение первого порядка du  f (u,t ) , dt для которого уравнение f (u,t )  0 имеет единственное решение u(0) t Рис. 10.3. К определению жесткости дифференциального уравнения u (t )  G (t ) . Любая интегральная кривая такого дифференциального уравнения характеризуется двумя участками с существенно различным поведением решения (см. рис. 10.3), причем первый участок значительно меньше второго. Первый участок с быстрым изменением функции u (t ) отраu жает стремление интегральной кривой к графику функции u (t )  G(t ) и называется пограничным слоем. На втором участке интегральная кривая практически совпадает с графиком G (t ) . Однако даже при небольшом отклонении от графика G (t ) в любой его точке производная du dt резко возрастает по сравнению с производной dG dt . Именно по этой причине малый шаг интегрирования, используемый при воспроизведении быстропротекающего участка, не может быть существенно увеличен вне пограничного слоя в явных методах численного интегрирования. Рассмотрим примеры жестких систем. 1. Простейшее дифференциальное уравнение du  u,   0 (10.1) dt 43 может быть жестким (рис. 10.4) , если интервал наблюдения [0, T ] знаu T t Рис. 13.4. Семейство решений уравнения (13.1) Рис. 10.4. Семейство решений уравнения (10.1) чительно превосходит величину   1 /  . Решением такого уравнения является функция u (t )  u (0)e t . Пунктирные линии на рис. 10.4 соответствуют различным значениям u (0) . 2. Для иллюстрации явления жесткости дифференциальных уравнений может быть полезной также любая система линейных ОДУ, матрица которой характеризуется большим разбросом собственных значений, например система Здесь собственные числа матрицы есть Решение этой системы (рис. 10.5) du1  1000u1  999u2 , dt du 2  u1  2u2 . dt 1  1001, 2  1 . (10.2) u1 (t )  0.999[u1 (0)  u2 (0)]e 1001t  [0.001u1 (0)  0.999u2 (0)]e t , u2 (t )  0.001[u1 (0)  u2 (0)]e 1001t  [0.001u1 (0)  0.999u2 (0)]e t содержит быструю составляющую как для u1 (t ) , так и для u2 (t ) , хотя по амплитуде быстрая составляющая функции u1 (t ) значительно превосходит аналогичную компоненту функции u2 (t ) (ее на рис. 10.5 в таком масштабе отобразить не удалось). u1,u2 u2(0) u1(t) u1(0) u2(t) t Рис. 10.5. Решение жесткой системы ОДУ (10.2) Завершая краткую характеристику свойства жесткости дифференциальных уравнений, отметим, что при решении научных задач жесткость уравнений является скорее правилом, чем исключением. По этой причине для таких задач разработаны специальные методы численного интегрирования. 44 Библиографический список рекомендуемых источников Основная литература 1. Ильин, В.П. Численные методы решения задач строительной механики : учеб.пособие для вузов / В.П.Ильин, В.В.Карпов, А.М.Масленников .— 2-е изд., доп.и перераб. — М.;СПб., 2005 .— 425 с. 2. Бахвалов Н.С. Численные методы.- М.: Наука, 2006. – 631 с. 3. Вержбицкий, В.М. Основы численных методов : учебник для вузов / В.М.Вержбицкий .— М. : Высш.шк., 2005 .— 840с. 4. Протасов И.Д. Лекции по вычислительной математике. - М.: Гелиос АРВ, 2009. – 211 с. 5. Бахвалов, Н.С. Московский государственный университет им. М.В. Ломоносова. Численные методы : учеб. пособие для вузов / Н.С. Бахвалов, Н.П. Жидков, Г.М. Кобельков .— 5-е изд. — М. : БИНОМ. Лаборатория Знаний, 2007 .— 636 с. Дополнительная литература 1. Кетков, Ю.Л. MATLAB 7:программирование,численные методы / Ю.Л.Кетков, А.Ю.Кетков, М.М.Шульц .— СПб. : БХВ-Петербург, 2005 .— 752 с. 2. Дьяконов, В. Maple 6 : Учеб.курс / В.Дьяконов .— СПб. и др. : Питер, 2001 . — 608 с. 3. Землянский, А.А. Практикум по информатике : учеб.пособие для вузов / А.А.Землянский, Г.А.Кретова,Ю.Р. Стратонович, Е.А.Яшкова;под ред.А.А.Землянский .— М. : КолосС, 2004 .— 384с. — (Учебники и учебные пособия для высших учебных заведений) .— Библиогр.в конце кн. — ISBN 5-9532-0046-3 /в пер./ : 179.41. 4. Микшина, В.С. Лабораторный практикум по информатике : Учеб.пособие для вузов / В.С.Микшина,Г.А.Еремеева,К.И.Бушмелева и др.;Под ред.В.А.Острейковского .— М. : Высш.шк., 2003 .— 376с. : ил. — Библиогр.в конце кн. — ISBN 5-06-004273-1 /в пер./ : 84.00. 45