Вычислительная математика лекции. Конспект лекций. Санкт петербург 2011 2 Оглавление
 Скачать 3.86 Mb.
|
|
8.8. Примеры использования квадратурных формул. Вычислить 3 1 2 dx x . Точное значение ln(5) = 1.6094379… 1) Используем составную формулу трапеций для n = 4 (число узлов равно 5). 64 h = (3- (-1) )/4 = 1 1 0 1 1 1 0 1 ( ) 2 2 2 1 1 / 5 1 1 1 1* 1.68333... 2 2 2 3 4 i i i n n n n i i h I f f f f I h f f f f h f Погрешность реальная |R| ист = 0.074. Число верных значащих цифр 2. Оценка погрешности 2 ' ' 2 2 ( ) ( ) max | ( ) | max | ( ) | 12 12 сост тр b a b a R f h f h M 2 2 3 2 1 2 ''( ) , M | ''( 1) | 2, |R| * 4 *1 * 2 (2 ) 12 3 f x f x Гарантирована только одна верная значащая цифра. Рассчитаем оценку h для заданной погрешности δ. 2 12 ( ) h b a M Если δ = 10 -2 , то h 0.122, а n 4 32.78. 0.122 b a h При n = 33, I = 1.61011… Вывод: при уменьшении h оценка погрешности становится более точной. Оценим погрешность, используя правило Рунге. Для этого x x 0 -1 x 1 0 x 2 1 x 3 2 x 4 3 f(x) 1 1/2 1/3 1/4 1/5 65 сформируем таблицу с шагом h/2=1/2/ x x 0 -1 x 1 -1/2 x 2 0 x 3 1/2 x 4 1 x 5 1.5 x 6 2 x 7 2.5 x 8 3 f(x) 1 2/3 ½ 2/5 1/3 2/7 ¼ 2/9 1/5 I=0.5((1+1/5)/2+2/3+1/2+1/3+2/7+1/4+2/9)=1.628968. Согласно правилу Рунге погрешность ΔI h/2 =1/3(1.6833-1.628968)=0.0195303. Относительно точного значения погрешность равна 1.628968-1.6094379=0.01446541 2) Используем составную формулу Симпсона. n = 4 четное. Воспользуемся формулой без половинных узлов. m =2. 1 0 2 1 2 1 1 1 1 1 1 1 ( 4 2 ) 1 4 2 1.62222... 3 3 5 2 4 3 m m n i i i i h I f f f f 4 (4) 4 4 ( ) ( ) max | ( ) | max | ( ) | 180 180 сост симп b a b a R f h f h M (4) 4 5 24 ( ) , M 24. (2 ) f x x |R| 4 24 8 * 4 * (1) 0.5333... 180 15 Оценка превышает соответствующую для формулы трапеций, потому что h = 1. 66 При заданном δ = 10 -2 2 4 180 180 *10 4 = = 0.37, n 10.8 ( ) 4 * 24 0.37 h b a M 8.9. Применение формул Ньютона - Котеса. 0 1 ( ) n н к k k k I b a H f N . Здесь n степень интерполяционного полинома, заменяющего подынтегральную функцию. Элементы H k , N рассчитываются в зависимости от n по формулам, полученным при анализе ИКФ. Значения этих элементов сведены в таблицу. m=n-1 H 0 H 1 H 2 H 3 H 4 H 5 H 6 N R m (f) 1 1 1 2 3 (2) ( ) ( ) 12 b a f 2 1 4 1 6 5 (4) ( ) ( ) 2880 b a f 3 1 3 3 1 8 5 (4) ( ) ( ) 6480 b a f 4 7 32 12 32 7 90 7 (6) ( ) ( ) 1935360 b a f 5 19 75 50 50 75 19 288 7 (6) 11( ) ( ) 87800000 b a f 6 41 216 27 272 27 216 41 840 Вычислим 3 1 2 dx x , используя 7-ординатную формулу Ньютона – Котеса. 4 2 6 6 3 b a h . Вычисления сводим в таблицу. 67 k x k f k H k H k f k 0 -1 1 41 41 1 -1/3 3/5 216 129.6 2 1/3 3/7 27 11.571 3 1 1/3 272 90.666 4 5/3 3/11 27 7.3636 5 7/3 3/13 216 49.846 6 3 1/5 41 8.2 338.2466 0 1 1 ( ) 4 338.2466 1.61704216. 840 n н к k k k I b a H f N Погрешность относительно точного значения ΔI=1.27*10 -3 Вычислим 1 4 2 1 ( 25 45 8) 4 x x dx с помощью 5- ординатной формулы Ньютона – Котеса. Эта формула должна быть точна, если подынтегральная функция есть полином степени не выше 4. Проверим это. h=(b-a)/n = 2/4=1/2.Вычисления сводим в таблицу как в предыдущем примере. k x k f k H k H k f k 0 -1 12 7 84 1 -1/2 1.6875 32 54 2 0 -8 12 -96 3 1/2 1.6875 32 54 4 1 12 7 84 Σ= 180 68 I=2*180/90 = 4. 8.10. Применение формул Гаусса. ИКФ Гаусса имеет вид 0 ( ) n G k k k I A f x , где n степень заменяющего интерполяционного полинома. Значения A k и x k в зависимости от n сведены в таблицу. Узлы n и веса 0 1 2 3 t 0 A 0 0.0000 2.0000 1 0.577350 3 1.0000 - √0.6 5/9 -0.861136 0.3478448 t 1 A 1 1 0.577350 3 1.0000 0.0000 8/9 -0.3399810 0.6521451543 t 2 A 2 √0.6 5/9 0.3399810 0.6521451541 t 3 A 3 0.861136 0.3478448 Осуществляя замену переменных, переходят к стандартному отрезку [-1,1]: 1 1 ( ) 2 2 2 b a b a a b b a f x dx f t dt 69 В нашем варианте 3 1 1 1 1 1 2 2 2 (1 2 ) dx dt x t Вычислим предыдущий интеграл с помощью трехординатной формулы Гаусса. 5 1 8 1 5 1 2 1.6026935... 9 9 2 (1 2*0) 9 2 (1 2 0.6) 2 (1 2 0.6 I Погрешность относительно точного значения ΔI=6.7*10 -3 Для четырехординатной формулы получаем I=1.60843…, для пятиординатной I=1.609288… 8.11. Замечания о точности квадратурных формул. Погрешность для составных квадратурных формул R=O(h k ), h=(b-a)/n. Чем больше k и чем меньше h, тем более точен результат. Отсюда не следует, что в конкретных случаях более грубая формула дает худший результат. Пример: f(x)=-8+45x 2 -25x 4 ; I= 1 1 ( ) 4 f x dx При h=1 a) составная формула трапеций I сост.тр =1/2f(-1)+f(0)+1/2f(1)=6- 8+6=4; б) составная формула Симпсона I сост.симп. =1/3(f(-1)+4f(0)+f(1))= =1/3(12-32+12)=8/3. Однако, если использовать формулу трапеций при шаге h/2=1/2, получим 70 I тр.сост. =1/2f(-1)+f(-1/2)+f(0)+f(1/2)+1/2f(1)= =6+1.6875-8+1.6875+6=7.375. Таким образом, полученное ранее точное значение есть результат случайного совпадения. Дальнейшее уменьшение шага будет уменьшать погрешность в обоих случаях, причем скорость уменьшения погрешности будет выше при использовании метода Симпсона. Численное интегрирование при фиксированном шаге может привести к неприемлемым погрешностям (особенно для осциллирующих функций). Так применение пятиординатного метода Ньютона – Котеса для функции, график которой изображен на нижнем рисунке, дает полностью искаженный результат. 0 x1 x2 x3 x4 x5 0 0 0 Использование адаптивных процедур в большинстве случаев позволяет избежать неадекватных решений. 9. Дополнительные элементы линейной алгебры и теории матриц 71 Теория матриц определяет основные классы математических моделей линейных непрерывных и дискретных систем управления. Излагаемый материал предполагает знакомство с курсом линейной алгебры, с арифметическими операциями над матрицами и понятием определителя матрицы. 9.1.Собственные числа и векторы. Преобразование y=Ax (x R n , A R n×n ) приводит к тому, что образ y испытывает вращение и изменение длины относительно прообраза x. Однако существуют особые ненулевые вектора e, называемые собственными, которые сохраняют направление и изменяют только длину, т. е. выполняется соотношение Ae = λe. Здесь λ скаляр, называемый собственным числом матрицы A. Собственные векторы являются решением системы линейных однородных уравнений (A-λE)e = 0. Известно, что эта система будет иметь ненулевые решения, если rank (A-λE) < n. В этом случае должно быть det (A-λE) = 0. Решив алгебраическое уравнение n-ой степени det (A-λE) = 0, найдем n значений собственных чисел λ (с учетом их кратности). Совокупность собственных значений называют спектром матрицы. Решая систему однородных линейных уравнений (A-λE)e = 0, для попарно различных λ i , вычисляют набор собственных векторов. Одному значению λ i может соответствовать один или несколько собственных векторов. Число собственных векторов для каждого λ i называют геометрической кратностью λ i в отличие от его алгебраической кратности. Доказано, что геометрическая кратность λ i не превышает его алгебраическую кратность. Таким образом, мы имеем набор собственных чисел (спектр матрицы) и соответствующий набор собственных векторов. 72 Доказано, что собственные векторы из этого набора линейно независимы. Обозначим α(λ) и γ(λ) алгебраическую и геометрическую кратности λ. Если для каждого λ α(λ) = γ(λ), то имеем дело с матрицей простой структуры. В противном случае матрица называется дефектной. Количество собственных векторов в наборе для матрицы простой структуры равно размеру матрицы. Иными словами для матрицы простой структуры количество собственных векторов равно количеству собственных значений с учетом их кратности. 9.2.Решение системы линейных однородных уравнений с вырожденной матрицей. Ранг матрицы. Пусть A R m×n . Произвольным образом выбираем "к" строк и "к" столбцов, где к min(m,n). Элементы, стоящие на пересечении этих строк и столбцов, образуют квадратную матрицу размером "к". Определитель этой матрицы называют минором к-ого порядка матрицы A. Определение. Рангом r матрицы называют максимальный порядок минора матрицы, отличного от нуля. Ранг нулевой матрицы равен нулю. Разность (min(m,n) – r) – дефект матрицы. Если дефект равен нулю, A m×n есть матрица полного ранга. Квадратную матрицу в этом случае называют невырожденной. Теорема. Если n число неизвестных однородной системы уравнений и r ранг её матрицы, то размерность пространства решений равна k = n –r (число ненулевых линейно независимых решений ). Если матрица полного ранга, то решение тривиально. 73 Далее полагаем, что A R n×n матрица квадратная. Алгоритм решения. Убеждаемся, что матрица А вырождена. Определяем ранг А, т. е. находим минор наибольшего размера, отличный от нуля. Из матрицы А выделяем подматрицу В r×n , состоящую из строк, выбранных при определении минора. Пусть при определении минора мы выбрали строки с номерами i 1 , i 2 ,…, i r и столбцы с номерами j 1 , j 2 ,…, j r . Будем решать систему Вх =0 следующим образом. Компоненты х с номерами, отличными от j 1 , j 2 ,…, j r (n –r элементов), будем считать произвольно заданными постоянными с 1 , с 2 ,…, с n-r . Тогда оставшиеся r неизвестных однозначно определятся из решения системы с квадратной невырожденной матрицей размера r. Число линейно независимых решений определяется количеством линейно независимых (n –r) мерных векторов (с 1 , с 2 ,…, с n-r ) T и равно n – r. Набор таких векторов удобно сформировать следующим образом: с (i) = 1 { } n r k k (i=1,…,n-r), 0 i 1 i=k k если k если Другой вариант. Воспользуемся методом исключения Гаусса. Метод включает два этапа. Сначала, используя последовательность эквивалентных преобразований, приводят матрицу системы к верхнему треугольному виду (прямой ход). На втором этапе решают полученную систему снизу вверх, начиная с последнего уравнения (обратный ход). На этапе прямого ходапоследовательно обнуляют элементы столбцов, расположенные ниже главной диагонали, 74 начиная с первого столбца и кончая(n-1) - ым. Для этого при обработке очередного i-ого столбца выделяют ненулевой элемент a ii , стоящий на главной диагонали матрицы, называемый ведущим элементом. Для обнуления каждого из элементов a ki , где k = i+1,…,n, из уравнения с номером k вычитают уравнение с номером i, умноженное на коэффициент a ki /a ii Если оказывается a ii = 0, среди элементов i-ого столбца ниже главной диагонали находят ненулевой элемент и переставляют соответствующее уравнение с i-ым уравнением. После обнуления элементов i-ого столбца переходят к следующему (i+1)-ому столбцу. В качестве очередного ведущего элемента пытаются использовать элемент a i+1 i+1 Если все анализируемые элементы i-ого столбца ниже главной диагонали оказались нулевыми (последнее возможно только при вырожденной матрице), говорят, что i-ый столбец не имеет ведущего элемента. Тогда в качестве нового ведущего элемента в i-ой строке пытаются использовать элемент a i i+1 в соседнем (i+1)-ом столбце и производят обнуление элементов (i+1) – ого столбца. После обработки всех (n-1) столбцов матрица будет иметь следующий вид. 1.Ненулевые строки расположены выше нулевых ( в противном случае нужно было произвести перестановку строк) и каждый ведущий элемент является первым ненулевым элементом всвоей строке. 2.В каждом столбце все элементы, расположенные ниже ведущего, равны нулю. 3.Каждый ведущий элемент лежит правее ведущих элементов предыдущих строк. 75 Пусть полученная матрица имеет r ведущих элементов, т. е. последние (n-r) строк - нулевые. Тогда ранг матрицы равен r, r переменных будут базисными, а (n-r) переменных свободными. Причем они соответствуют столбцам преобразованной матрицы с ведущими элементами и без них соответственно. Компоненты вектора переменных, отнесенные к свободным переменным, задаются набором произвольных констант с 1 , с 2 ,…, с n-r . Базисные переменные вычисляются через свободные переменные с помощью обратной подстановки (обратный ход). Число линейно независимых решений определяется количеством линейно независимых (n –r) мерных векторов (с 1 , с 2 ,…, с n-r ) T и равно n – r. Набор таких векторов удобно сформировать следующим образом: с (i) = 1 { } n r k k (i=1,…,n-r), 0 i 1 i=k k если k если . Пример. Решить однородную систему Ax=0, 1 1 5 1 1 1 2 3 3 1 8 2 1 3 9 7 A . Осуществляем процедуру Гауссова исключения A ⟶A 1 ⟶A 2 ⟶U. 1 1 1 5 1 0 2 7 4 0 2 7 5 0 4 14 8 A 2 1 1 5 1 0 2 7 4 0 0 0 1 0 0 0 0 A Обработка третьего столбца не потребовалась, A 2 = U. r(A) =r(U)=3, n –r = 1. Существует одно ненулевое линейно независимое решение. Третий столбец матрицы U оказался без ведущего элемента. Переменная x 3 =c свободная переменная. 76 Остальные переменные находятся обратной подстановкой. Из третьего уравнения следует x 4 = 0, из второго 2x 2 -7c + 4x 4 = 0, т. е. x 2 = (7/2)c. Из первого уравнения имеем x 1 – x 2 + 5x 3 – x 4 = 0, т. е. x 1 = -(3/2)c. 3 3 2 2 7 7 2 2 1 0 0 c c x c c 9.3.Подобное преобразование . Говорят, что квадратные матрицы A и B подобны, если существует невырожденная матрица P такая, что B = P -1 AP. Подобные матрицы имеют одинаковый набор собственных чисел. Действительно det(P -1 AP – λE) = det(P -1 (A – λE)P) = =detP -1 det(A – λE)detP= det(A – -λE). Преобразование подобия возникает естественным путем как результат замены переменных (или перехода к новому базису) в пространстве n мерных векторов. Действительно, пусть модель процесса (устройства) имеет вид y=Ax. Произведем замену переменных x=Px1, y=Py1. Тогда уравнение y=Ax примет вид y1= P -1 APx1.Это означает, что в новых переменных то же самое преобразование осуществляется уже не матрицей A, а матрицей P -1 AP, подобной A. Подобное преобразование применяют для перехода к новой матрице более простой формы, пригодной для решения задачи. Например, если удастся путем последовательности подобных преобразований привести исходную матрицу к треугольному виду, будут вычислены собственные значения, стоящие на диагонали этой матрицы. |