Курсовая работа по численным методам для ii курса методические указания СанктПетербург 2009
 Скачать 317.74 Kb.
|
|
САНКТ-ПЕТЕРБУРГСКИЙ ГОСУДАРСТВЕННЫЙ УНИВЕРСИТЕТ Факультет прикладной математики – процессов управления А. П. ИВАНОВ, И. В. ОЛЕМСКОЙ КУРСОВАЯ РАБОТА ПО ЧИСЛЕННЫМ МЕТОДАМ ДЛЯ II КУРСА Методические указания Санкт-Петербург 2009 Перейти к оглавлению на странице: 26 КУРСОВАЯ РАБОТА ПО ЧИСЛЕННЫМ МЕТОДАМ 1. Методом Рунге-Кутты IV порядка проинтегрировать диффе- ренциальное уравнение dy dx + a · y = ϕ(x), ϕ(x) = sin kx, k = aπ 4 , (1) на отрезке x ∈ [0; 5] , приняв π = 3.14159265 , с начальным условием y(0) = 0 . Вывести значения интегрируемой функ- ции в точках x i = 0.5 · i, i = 0, 10. 2. Используя представление искомого решения в виде y(x) = Z x 0 e a(t−x) ϕ(t) dt, (2) получить значения y(x) с точностью ε в заданных точках x i = 0.5 · i , i = 0, 10 , применив для вычисления интеграла составную квадратурную формулу Симпсона. 3. Найти с заданной точностью ε ближайший к нулю корень уравнения y(x) = 0 . 4. Взяв значения z i = y(i) + (−1) i ε i , i = 0, 10 , построить полином наилучшего среднеквадратичного приближения пя- той степени ¯ Q 5 (x) , то есть найти параметры {a i } полинома Q 5 (x) = a 0 x 5 + a 1 x 4 + a 2 x 3 + a 3 x 2 + a 4 x + a 5 из условия 10 X i=0 ( ¯ Q 5 (x i ) − z i ) 2 = min {a i } 10 X i=0 (Q 5 (x i ) − z i ) 2 (3) Решить эту же задачу, заменив {z i } точными значениями функции y(x) в точках x i = i, i = 0, 5 (см. формулу (1.7) ). В заданиях 1, 2, 4 построить графики соответствующих функций. 2 Перейти к оглавлению на странице: 26 ГЛАВА 1. ЧИСЛЕННЫЕ МЕТОДЫ РЕШЕНИЯ ЗА- ДАЧИ КОШИ ДЛЯ ОБЫКНОВЕННЫХ ДИФФЕРЕНЦИАЛЬНЫХ УРАВНЕНИЙ §1. Постановка задачи Пусть в области D = {a 6 x 6 b, |y i − y i 0 | 6 b i } ∈ R n+1 определена функция f ≡ f (x, y 1 , . . . , y n ), (x, y) ∈ D. Необходимо найти решение dy dx = f (x, y) (1.1) удовлетворяющее начальному условию y(x 0 ) = y 0 (1.2) §2. Одношаговые методы Одношаговые методы – это методы, которые последователь- но дают приближения y m к значениям точного решения y(x m ) в каждом узле x m сетки на основе известного приближения y m−1 к решению в точке x m−1 В общем виде их можно представить : y m+1 = F (f, x m+1 , x m , y m+1 , y m ). (2.1) Надо отметить, что одношаговые методы вида (2.1) в правой ча- сти содержат искомое значение y m+1 . Это позволяет ввести еще один параметр классифицирующий численные методы. Наличие в правой части представления метода в виде (2.1) зависимости от искомой функции y m+1 делает его неявным. Отсутствие этой за- висимости – явным, вида: y m+1 = F (f, x m+1 , x m , y m ). (2.2) 3 Перейти к оглавлению на странице: 26 Общая схема явного m – этапного одношагового метода Рунге- Кутты имеет вид y(x + h) ≈ y(x) + m X i=1 p i k i (h), (2.3) k i (h) = hf (x + α i h, y(x) + i−1 X j=1 β ij k j ), i = 1, m. (2.4) Он идеально приспособлен для практического расчета: во-первых, не требует вычисления дополнительных начальных значений; во- вторых, позволяет легко менять шаг интегрирования. Функция Ψ(h) = % x+h = y(x + h) − y(x) − m X i=1 p i k i (h) ≈ O(h s+1 ), (2.5) называется методической (локальной) погрешностью одно- шагового метода (методической погрешностью на шаге), а s – порядком метода. Выбор постоянных (параметров метода) p i , α i , β ij произво- дится так, чтобы разложение методической погрешности (2.5) по степеням h начиналось со степени s + 1 при произвольной функ- ции f (x, y). Расчётная схема четырёхэтапного метода Рунге- Кутты четвёртого порядка. (Правило одной шестой). y(x + h) ≈ y(x) + 1 6 [k 1 (h) + 2k 2 (h) + 2k 3 (h) + k 4 (h)]. (2.6) где k 1 (h) = hf (x, y(x)), k 2 (k) = hf (x + 1 2 h, y(x) + 1 2 k 1 (h)), k 3 (k) = hf (x + 1 2 h, y(x) + 1 2 k 2 (h)), k 4 (k) = hf (x + h, y(x) + k 3 (h)). Погрешность расчётной схемы (2.6) , как и для всех формул Рунге- Кутты четвёртого порядка точности, представляется в следующей 4 Перейти к оглавлению на странице: 26 форме: % x+h = Ψ (5) (0) 120 h 5 + o(h 5 ). (2.7) Классификация погрешностей. Требуется найти функцию y(x), которая является решением диф- ференциального уравнения dy dx = f (x, y(x)), x ∈ [a, b], (2.8) и принимает в точке x 0 = a некоторое определённое значение y(x 0 ) = y 0 (2.9) • В случае, если начальное условие y(x 0 ) вычислено (задано) с погрешностью, то вместо задачи (2.8) , (2.9) , решается задача: dy 0 dx = f (x, y 0 (x)), x ∈ [a, b], (2.10) y 0 (x 0 ) = ¯ y 0 (2.11) с изменёнными начальными условиями y 0 − ¯ y 0 = R 0 6= 0. (2.12) Решение задачи (2.10) , (2.11) зависит от ¯ y 0 и не совпадает с решением y(x) исходной задачи (2.8) , (2.9) Разность ξ n = y(x n ) − y 0 (x n ) называется неустранимой по- грешностью решения y 0 (x) . • Разность между значением решения y 0 (x n ) задачи (2.10) , (2.11) и его приближённым значением y n , полученным по формуле y n = F (f ; h; x n−1 ; y n−1 ) (2.13) одношагового метода, ε n = y 0 (x n ) − y n (2.14) 5 Перейти к оглавлению на странице: 26 называется погрешностью метода, или глобальной погрешно- стью метода. Но вследствие ошибок округления и приближённого вычисле- ния правой части f (x, y) дифференциального уравнения вы- числения значений y n по формуле (2.13) выполняются неточ- но. Фактически найденные значения ˆ y n удовлетворяют не соотношению (2.13) , а условию F (f ; h; x n−1 ; ˆ y n−1 ) − ˆ y n = δ n (2.15) Невязка δ n называется погрешностью округления на n -ом шаге. • Разность между точным решением y(x n ) задачи (2.8) , (2.9) и приближённым фактически найденным значением ˆ y n R n = y(x n ) − ˆ y n (2.16) называется полной погрешностью приближённого ре- шения • Величина η n = y n − ˆ y n называется вычислительной по- грешностью. • Из соотношений (2) , (2.13) , (2.16) , (2) следует, что R n = ξ n + ε n + η n (2.17) т.е полная погрешность приближённого решения рав- на сумме неустранимой погрешности, погрешности метода и вычислительной погрешности. Для полной погрешности справедливо представление R n = n X j=1 (% j + δ j ) · e R xn xj f y (ψ,¯ y j (ψ)) dψ + R 0 · e R xn x0 f y (ψ,¯ y 0 (ψ)) dψ (2.18) Здесь y j (x) – интегральная кривая, проходящая через точку (x j , ˆ y j ). 6 Перейти к оглавлению на странице: 26 Из (2.18) следует, что полная погрешность приближённого ре- шения задачи Коши (2.8) , (2.9) в точке x n равна сумме локальных погрешностей на каждом шаге, взятых с коэффициентами e R xn xj f y (ψ,¯ y j (ψ)) dψ (2.19) Характер отклонения приближённого решения от точного, т.е. эво- люция полной погрешности, зависит от поведения интегральных кривых уравнения (2.8) • Если ∂f ∂y = 0, т.е. правая часть уравнения не зависит от y ,то полная погрешность равна сумме локальных погрешностей: R n = R 0 + n X j=1 (% j + δ j ) (2.20) • Если ∂f ∂y > 0, т.е. интегральные кривые расходятся, влияние локальных погрешностей, полученных на предыдущих шагах, возрастает и полная погрешность больше суммы локальных погрешностей. • Если ∂f ∂y < 0, т.е. интегральные кривые сближаются, влияние локальных погрешностей, полученных на предыдущих шагах, ослабевает и полная погрешность меньше суммы локальных погрешностей. Такое поведение характерно как для полной погрешности, так и для отдельных её частей. 7 Перейти к оглавлению на странице: 26 Практическая реализация явных одношаговых мето- дов типа Рунге-Кутты решения задачи Коши. Полагаем, что в нашем вычислительном процессе погрешно- стью округления и неустранимой погрешностью можно пренебречь. Оценки погрешностей необходимы, с одной стороны, чтобы обеспечить длину шага h, достаточно малую для достижения тре- буемой точности вычисляемых результатов, а с другой стороны – чтобы гарантировать достаточно большую длину шага во избежа- ние бесполезной вычислительной работы. Метод Рунге оценки полной погрешности. Полагаем, что в точке x n по m -этапному методу ( 2.3 ) s -ого порядка точности с постоянным шагом h вычислено приближён- ное решение ¯ y n исходной задачи Коши. Считаем, что для полной погрешности метода Рунге-Кутты s− о порядка справедливо пред- ставление: y(x n ) − ¯ y n ≈ z(x n )h s , (2.21) где z(x n )h s − главный член аcимптотического разложения полной погрешности. Используя ту же расчётную формулу с шагом h 2 , вычислим в той же точке x n другое значение решения ˜ y n . Для этого потребуется в два раза больше шагов и погрешность прибли- жённого решения в этом случае может быть представлена y(x n ) − ˜ y n ≈ z(x n ) h 2 s (2.22) Исключим из (2.21) , (2.22) точное значение решения y(x n ) И пра- вило Рунге оценки глобальной погрешности может быть записано следующим образом: ¯ R n = y(x n ) − ¯ y n ≈ (˜ y n − ¯ y n ) (1 − 2 −s ) (2.23) ˜ R n = y(x n ) − ˜ y n ≈ (˜ y n − ¯ y n ) (2 s − 1) (2.24) В качестве решения в точке x n примем значение ˜ y n как более точное по сравнению с ¯ y n . Для него имеем оценку погрешности 8 Перейти к оглавлению на странице: 26 (2.24) . Эта величина может быть как больше, так и меньше некото- рого значения ε, являющегося наперёд заданной допустимой погрешностью. Если | ˜ R n | 6 ε, то заданная точность приближённого решения достигается. Полученное приближённое решение можно уточнить, прибавив к нему величину главного члена погрешности, т.е. поло- жив y(x n ) ≈ y n = ˜ y n + ˜ R n При этом y(x n ) − y n = O(h s+1 ). 9 Перейти к оглавлению на странице: 26 ГЛАВА 2. ВЫЧИСЛЕНИЕ ОПРЕДЕЛЕННОГО ИНТЕ- ГРАЛА §1. Квадратурные формулы Будем рассматривать задачу о вычислении однократного ин- теграла J (F ) = Z b a F (x) dx (1.1) с помощью конечного числа значений интегрируемой функции. Интервал интегрирования [a, b] есть любой конечный (или бесконечный) отрезок числовой оси. Подынтегральная функция F (x) – любая интегрируемая в смысле Римана функция. Каждое правило основано на замене интегрируемой функции на какую-либо элементарную функцию – алгебраический много- член, рациональную функцию, тригонометрический многочлен и др. Чтобы такая замена давала хорошую точность, требуется, что- бы заменяемая функция F (x) обладала высоким порядком глад- кости. Если F (x) имеет какие-нибудь особенности, мы заинтересо- ваны в выделении их. Выделение делается обычно при помощи разложения F (x) на два сомножителя F (x) = p(x) · f (x), где • p(x) имеет особенности того же типа, что и F (x), и назы- вается весовой функцией ; • f (x) – гладкая функция. Такое разложение приводит (1.1) к виду: J (F ) = Z b a p(x)f (x) dx. (1.2) Считая вес p(x) фиксированным, а f (x) любой гладкой функцией, будем строить правило интегрирования, рассчитанное на функции, имеющие одинаковые, заранее известные особенности. С учётом 10 Перейти к оглавлению на странице: 26 всего вышесказанного будем строить правила вычисления опреде- лённого интеграла следующего вида: J (F ) = Z b a p(x)f (x) dx ≈ n X k=1 A k f (x k ) = S n , x k ∈ [a, b], (1.3) где • f (x) ∈ Φ(a, b) – некоторый класс функций f (x) , определён- ных на [a, b] ; • p(x) – весовая функция – некоторая фиксированная неот- рицательная на [a, b] функция, для которой Z b a p(x) dx > 0, (1.4) и для ∀f (x) ∈ R существует Z b a p(x)|f (x)| dx. (1.5) Определение 1.1. Равенство (1.3) называют формулой механи- ческих квадратур или просто квадратурной формулой, а S n = P n k=1 A k f (x k ) – квадратурной суммой; A k – квадратурными коэф- фициентами формулы (1.3) ; x k – узлами квадратурной формулы (1.3) Будем предполагать, что узлы квадратурной формулы (1.3) упорядочены по возрастанию a = x 0 < x 1 < . . . < x n = b , n, x k , A k , k = 1, . . . , n; – являются параметрами квадратур- ной формулы (1.3) и их следует выбирать так, чтобы достигнуть “возможно лучшего” результата интегрирования для всех функций избранного класса. Определение 1.2. Величина R n (f ) = Z b a p(x)f (x) dx − n X k=1 A k f (x k ) = J (f ) − S n (f ). (1.6) называется методической погрешностью (остаточным членом) квадратурной формулы (1.3) 11 Перейти к оглавлению на странице: 26 При зтом R n (f ) зависит от свойств функции f и от выбора квадратурной формулы, т. е. от узлов и коэффициентов. Определение 1.3. Квадратурная формула (1.3) имеет алгебраиче- скую степень точности m, если она верна для любых многочленов степени m и не верна для многочленов степени m + 1 . Это равносильно тому, что равенство Z b a p(x)x i dx = n X k=1 A k x i k , (1.7) выполняется для i = 0, 1, . . . , m и не выполняетсядля i = m + 1. Или, что то же самое, R n (x i ) = 0, i = 0, 1, . . . , m; R n (x m+1 ) 6= 0. (1.8) §2. Квадратурная формула Симпсона Трёхточечная квадратурная формула – формула Симпсона – имеет вид : Z b a f (x) dx ≈ (b − a) 6 f (a) + 4 · f ( a + b 2 ) + f (b) (2.1) Она является точной для всех многочленов третьей степени, т.е. ее алгебраическая степень точности равна трем. В классе C 4 (a,b) остаточный член квадратурной формулы (2.1) R 3 (f ) = Z b a f (4) (ξ(x)) 4! (x − a)(x − a + b 2 ) 2 (x − b) dx = = − 1 90 ( b − a 2 ) 5 f (4) (η), η ∈ [a, b]. (2.2) Составные квадратурные формулы. Основная идея мето- да заключается в том, что для повышения точности интегрирова- ния отрезок [a, b] делят на несколько частей. Применяют избран- ную квадратурную формулу к каждой отдельной части и резуль- таты складывают. 12 Перейти к оглавлению на странице: 26 Для большинства квадратурных формул методическая по- грешность R n (f ) зависит от величины отрезка интегрирования и может быть представлена в виде: R n (f ) = (b − a) k C(a, b), (2.3) где C(a, b) есть медленно меняющаяся функция от a, b и k. Эта зависимость показывает, что если уменьшить отрезок интегри- рования в p раз, то R n (f ) уменьшится в p k раз. Чем больше k, тем значительнее будет уменьшение. Разобъём исходный интервал интегрирования на частичные отрезки: [a, b] : a = z 0 < z 1 < z 2 < . . . < z k = b. В самом общем случае z i+1 − z i = h i > 0 и h i предполагаются различными. Теперь применим на каждом из частичных отрезков [z i , z i+1 ] квадратурное правило (в общем случае своё) для вычисления ин- теграла: J (i) (f ) = Z z i+1 z i p(x)f (x) dx = n i X j=1 A ij f (x ij ) + R n i (f ) = = S i n i (f ) + R n i (f ), i = 0, k − 1. (2.4) Просуммируем по i левые и правые части (2.4) : J (f ) = Z b a p(x)f (x) dx = k−1 X i=0 J (i) = k−1 X i=0 Z z i+1 z i p(x)f (x) dx ≈ ≈ k−1 X i=0 n i X j=1 A ij f (x ij ) = k−1 X i=0 S i n i (f ) = S N (f ), (2.5) N = k−1 X i=0 n i ; R N = J (f ) − S N (f ) = k−1 X i=0 R n i (f ), R n i (f ) = J i (f ) − S i n i (f ). В качестве квадратурных формул (2.4) могут быть использованы квадратурные формулы всех существующих типов. Причём как в 13 Перейти к оглавлению на странице: 26 смешанном составе – в составной квадратурной формуле использу- ется несколько типов, так и в однородном – в составной квадратур- ной формуле используются только квадратурные формулы одного типа. Тоже самое можно сказать и о числе узлов квадратурной фор- мулы используемом на каждом из частичных отрезков - их число может зависеть от номера отрезка. Такой подход увеличения точ- ности применим к простейшим формулам Ньютона-Котеса (здесь – Симпсона). Составная квадратурная формула Симпсона. Разобъём интервал интегрирования [a, b] на k = 2m частичных отрезков, z i = a + iH, i = 0, k H = b−a k и будем применять квадратурное правило на сдвоенном частичном отрезке [a + (2j − 2)H, a + (2j)H] j = 1, m. В качестве квадратурного правила используемого на каждом из сдвоенных частичных отрезков будет использоваться квадратурное правило Симпсона (2.1) : J (j) (f ) = Z z 2j z 2j−2 f (x) dx ≈ z 2j − z 2j−2 6 [f (z 2j−2 ) + 4f (z 2j−1 )+ +f (z 2j )] = H 3 [f (z 2j−2 + 4f (z 2j−1 ) + f (z 2j )] , j = 1, m. (2.6) Просуммировав по всем сдвоенным частичным отрезкам [a, a + 2H], [a + 2H, a + 4H], . . . выпишем составную квадратурную фор- мулу парабол Z b a f (x) dx = b − a 3k [f (z 0 ) + f (z k ) + 2(f (z 2 ) + f (z 4 ) + . . . +f (z k−2 )) + 4(f (z 1 ) + f (z 3 + . . . + f (z k−1 ))] + R k (f ).(2.7) Воспользуемся результатами оценки методической погрешности для квадратурной формулы парабол (2.2) : R k+1 (f ) = − H 4 180 Z b a f (4) (x) dx + o(H 4 ) = (2.8) = C 4 H 4 + o(H 4 ).; o(H 4 ) H 4 → 0, H → 0, где C 4 ≡ C 4 (f ) ≡ − 1 180 Z b a f (4) (x) dx. 14 Перейти к оглавлению на странице: 26 Практические способы оценки погрешности составных квадратурных формул. Оказывается, что подобно представле- нию (2.3) для погрешности составной квадратурной формулы (2.5) для широкого класса функций f (·) можно написать: R k (f ) = C m H m + o(H m ). (2.9) где m – натуральное число, C m = C m (f ) – некоторая констан- та, зависящая лишь от f (·) и типа формулы, но не зависящая от H, o(H 4 ) H 4 → 0, H → 0. Формула (2.9) дает ассимптотиче- ское представление (разложение) погрешности (2.5) по параметру H - шагу равномерного разбиения интервала интегрирования. Она справедлива, в частности, когда малые квадратурные формулы (2.1) (Симпсона, например,) имеют алгебраическую степень точно- сти m − 1, а функция f (·) – непрерывную (интегрируемую) на [a, b] производную f m (·). Практическая оценка константы C m проводится следующим образом. Фиксируем два значения шага разбиения – H 1 = H и H 2 = H/L, L > 1, и проводим расчёты на двух равномерных сет- ках с шагом H 1 и H 2 соответственно по исследуемой квадратур- ной формуле. Предполагаем, что для её методической погрешности верно аcимптотическое разложение (2.9) , т.е. R (H 1 ) k (f ) = J (f ) − S (H 1 ) k = C m H m + o(H m ) (2.10) и R (H 2 ) k (f ) = J (f ) − S (H 2 ) k = C m H L m + o H L m (2.11) Исключая из (2.10) , (2.11) неизвестное значение интеграла J, най- дём S (H 2 ) k − S (H 1 ) k = C m H m − H L m + o(H m ) (2.12) а отсюда C m = S (H 2 ) k − S (H 1 ) k H m (1 − L −m ) + o(H m ) H m (2.13) 15 Перейти к оглавлению на странице: 26 Подставляя это выражение для C m в равенства (2.10) , (2.11) и от- брасывая бесконечно малые o(H m ), найдём приближённые пред- ставления погрешностей R (H 1 ) k , R (H 2 ) k составной квадратурной формулы на двух равномерных сетках : R (H 1 ) k (f ) = J (f ) − S (H 1 ) k ≈ S (H 2 ) k − S (H 1 ) k 1 − L −m (2.14) R (H 2 ) k (f ) = J (f ) − S (H 2 ) k ≈ S (H 2 ) k − S (H 1 ) k L m − 1 (2.15) Этот способ оценивания методической погрешности называется правилом Рунге. Так при L = 2, т.е. шаг дробления умень- шается в два раза, погрешность составной квадратурной формулы R (H 2 ) k (f ) с шагом H 2 меньше погрешности R (H 1 ) k (f ) составной квадратурной формулы с шагом H 1 почти в 2 m раз. Используя правило Рунге можно решить вопрос о нахожде- нии оптимального шага разбиения H o интервала интегрирования, обеспечивающего (в первом приближении) вычисление интеграла с заданной точностью ε. ε = |R (H o ) k (f )| ≈ |C m H m o | = S (H 2 ) k − S (H 1 ) k H m (1 − L −m ) H m o , (2.16) а отсюда и получаем порядок оптимального шага разбиения H o ≈ H ε(1 − L −m ) S (H 2 ) k − S (H 1 ) k 1 m (2.17) 16 Перейти к оглавлению на странице: 26 ГЛАВА 3. МЕТОД НЬЮТОНА РЕШЕНИЯ УРАВНЕ- НИЙ §1. Основы метода Рассмотрим вопрос о решении скалярного уравнения f (x) = 0, x ∈ [α, β] (1.1) с использованием метода Ньютона в предположении, что решение ¯ x уравнения существует и ¯ x ∈ [α, β] . При условии, что функция f (·) дифференцируема на рассматриваемом интервале и при на- личии хорошего приближения x k к корню ¯ x функции можно ис- пользовать метод Ньютона, называемый также методом линеари- зации или методом касательных. Его расчётные формулы могут быть получены заменой исходного уравнения (1.1) в окрестности корня линейным уравнением (здесь x k – приближение корня ¯ x ) f (x k ) + f 0 (x k )(x − x k ) = 0. (1.2) Решение этого уравнения принимается за очередное приближение x k+1 к искомому корню ¯ x исходного уравнения (1.1) : x k+1 = x k − f (x k ) f 0 (x k ) (1.3) Метод Ньютона имеет простую геометрическую интерпретацию: график функции заменяется касательной к нему в точке (x k , f (x k )) и за очередное приближение x k+1 принимается абсцисса точки пе- ресечения этой касательной с осью Ox . Используя эту интерпре- тацию, легко получить расчётные формулы (1.3) метода Ньютона и вследствие этой интерпретации он именуется, как указано выше, методом касательных. Ясно, что сходимость последовательности {x k } к корню зави- сит от свойств функции f (·) и не всегда имеет место. Пусть выбра- но начальное приближение x 0 . Легко представить, что уже прибли- жение x 1 не попадает на исходный интервал определения функции и процесс итераций останавливается. Приведём полезную теорему, гарантирующую сходимость ме- тода Ньютона. 17 Перейти к оглавлению на странице: 26 Теорема 1.1. Если для функции f (·) ∈ C 2 [α, β] выполнены усло- вия • f (α) · f (β) < 0 , • f 0 (x) 6= 0 и f 00 (x) 6= 0 на [α, β] (и, следовательно, сохраняют определённые знаки при x ∈ [α, β]) , то исходя из начального приближения x 0 ∈ [α, β] , для которого f (x 0 ) · f 00 (x 0 ) > 0 , по формуле (1.3) можно вычислить единствен- ный корень ¯ x уравнения (1.1) с любой степенью точности. ♣ Замечание. Практическим критерием окончания вычислений является выполнение условия |x n+1 − x n | < ε , где ε – требуемая точность вычисления корня. ♣ Если же условия теоремы 1.1 не выполняются или проверка их затруднительна, то очередное “приближение” x k+1 может оказать- ся вне интервала, на котором расположен корень ¯ x . В этом случае x k+1 строится либо методом половинного деления либо методом хорд. В первом случае полагают x k+1 = a k + b k 2 , (1.4) во втором – x k+1 = a k − b k − a k f (b k ) − f (a k ) · f (a k ). (1.5) Здесь a k , b k – левый и правый конец интервала, которому принад- лежит корень ¯ x на предыдущем шаге. На начальном этапе полагаем a 0 = α, b 0 = β . Пусть для определённости f (a) < 0, f (b) > 0. Если x 1 ∈ [α, β] , то вычислив c = f (x 1 ) , полагаем a 1 = c , b 1 = b 0 при c < 0 , и a 1 = a 0 , b 1 = c при c > 0 и повторяем вычисления. Если же приближение x 1 6∈ [α, β] , то применяем формулы (1.4) либо (1.5) и поступаем как и выше: вычисляя c = f (x 1 ) , полагаем a 1 = c , b 1 = b 0 при c < 0 , и a 1 = a 0 , b 1 = c при c > 0 и применяем метод Ньютона. 18 Перейти к оглавлению на странице: 26 Комментарии к выполнению задания Формула (1.3) для функции, заданной уравнением (1) , запи- шется в виде x k+1 = x k − y(x k ) y 0 (x k ) = x k − y(x k ) ϕ(x k ) − ay(x k ) (1.6) Следовательно, для её использования достаточно проинтегриро- вать уравнение (1) на отрезке [0, x k ] , либо вычислить интеграл (2) в пределах от 0 до x k . Можно также использовать аналитическое представление для решения уравнения (1) , которое имеет вид y(x) = a · sin kx − k · cos kx + ke −ax a 2 + k 2 (1.7) Начальное приближение для использования формулы (1.6) может быть найдено, например, графическим методом, т.к. к моменту ре- шения этой задачи уравнение (1) уже проинтегрировано. ♣ 19 Перейти к оглавлению на странице: 26 ГЛАВА 4. ЗАДАЧА АППРОКСИМАЦИИ Для решения поставленной в пункте 4 задачи использу- ется метод наименьших квадратов. он состоит в следующем. §1. Линейная задача метода наименьших квадратов Пусть на отрезке [a, b] задана некоторая функция f (·) , причём её значения {y i = f (x i )} в узлах сетки {x 0 , x 1 , . . . , x n } известны с погрешностями {ε i } , то есть вместо набора значений {y i } имеем набор {˜ y i = y i +ε i } . (Далее под {y i } будем понимать заданные зна- чения функции, т.е. с погрешностями, вводя обозначение ˜ y i лишь в случае необходимости). Пусть, кроме того, на [a, b] определены функции ϕ j (·) ∈ Φ , j = 0, m . Введём в рассмотрение обобщённый полином P m (x) = m X j=0 a j ϕ j (x). Пусть a – вектор коэффициентов полинома P m (·) , y – вектор зна- чений функции f (·) и, наконец, ϕ(·) – вектор-функция, составлен- ная из значений {ϕ i (x)} : a = (a 0 , a 1 , . . . , a m ) T , y = (y 0 , y 1 , . . . , y n ) T , ϕ(x) = (ϕ 0 (x), ϕ 1 (x), . . . , ϕ m (x)) T , Введём также функции σ(a, y) = n X i=0 (P m (x i ) − y i ) 2 , δ(a, y) = r σ(a, y) n + 1 Функцию δ(a, y) назовём среднеквадратичным уклонением обоб- щённого полинома P m (·) от функции f (·) на системе узлов {x i } . Отказываясь теперь от условия P m (x i ) = y i , i = 0, n , поставим задачу так: найти обобщённый полином ¯ P m (·) = ¯ a T ϕ(·) , для кото- рого среднеквадратичное уклонение минимально : δ(¯ a, y) = min a δ(a, y). 20 Перейти к оглавлению на странице: 26 Поставленную здесь задачу называют линейной задачей метода наименьших квадратов или просто методом наименьших квадра- тов (МНК). Если искомый полином существует, то будем называть его многочленом наилучшего среднеквадратичного приближения. Отметим, что минимум функций σ(·) и δ(·) достигается на одном и том же векторе ¯ a , поэтому фактически будем вести ми- нимизацию функции σ(a, y) , а не δ(·) . Простейший подход к решению задачи – использование необ- ходимых условий в задаче поиска экстремума для дифференциру- емой функции: ∂σ(a, y) ∂a k a=¯ a = 0, k = 0, m. (1.1) В дополнение к уже введённым ранее обозначениям примем также следующее: Q = ϕ 0 (x 0 ), ϕ 1 (x 0 ), ϕ m (x 0 ), ϕ 0 (x 1 ), ϕ 1 (x 1 ), ϕ m (x 1 ), ϕ 0 (x n ), ϕ 1 (x n ), ϕ m (x n ) Кроме того будем считать, что под скалярным произведением двух векторов u, v ∈ R k понимается число (u, v) = u T v = u · v = k X i=1 u i v i , а под нормой векторов – евклидова норма kyk 2 = (y, y) . Тогда в новых обозначениях σ(a, y) = kQa − yk 2 = (Qa − y, Qa − y) = = (Qa, Qa) − 2(Qa, y) + (y, y) = = (a, Q T Qa) − 2(a, Q T y) + kyk 2 Для определения параметров искомого полинома в соответствии с формулой (1.1) имеем СЛАУ: Ha = b, где H = Q T Q, b = Q T y. (1.2) 21 Перейти к оглавлению на странице: 26 В предложенной в пункте 4 задаче её параметры соотносятся с ис- пользованными выше параметрами следующим образом: a = (a 0 , a 1 , . . . , a 5 ) T , y = z = (z 0 , z 1 , . . . , z 10 ) T , ϕ(x) = (1, x, . . . , x 5 ) T , σ(a, y) = σ(a, z) = 10 X i=0 (Q 5 (x i ) − z i ) 2 , δ(a, z) = r σ(a, z) 5 , Q = 1, x 0 , x 5 0 , 1, x 1 , x 5 1 , 1, x 5 , x 5 5 Остаётся для определения параметров полинома ¯ Q 5 (x) решить си- стему линейных алгебраических уравнений (1.2) . В предложенном варианте матрица H системы (1.2) неособая и задача имеет един- ственное решение. Выполняя вторую часть задания из п.4, Вы должны получить так называемый интерполяционный полином ˜ Q 5 (x) , обладающим свойством: ˜ Q 5 (i) = y(i) , i = 0, 5 . (Построив полином – убедиться!) §2. Методы Гаусса Метод Гаусса относится к точным методам решения систем линейных алгебраических уравнений (СЛАУ). Под точными мето- дами решения СЛАУ понимают методы, которые позволяют за ко- нечное число арифметических операций при отсутствии округле- ний получить точное решение СЛАУ при точно заданной правой части СЛАУ и и точно заданных элементах её матрицы. Перечис- лим здесь некоторые из них с указанием алгоритмов их использо- вания. Методы Гаусса (или методы исключения неизвестных) опти- мальны для решения СЛАУ общего вида по количеству арифмети- ческих операций, необходимых для нахождения решения этой си- стемы. Пусть СЛАУ записана в виде Ax = b. (2.1) Простой метод Гаусса состоит в следующем: в СЛАУ (2.1) , запи- санной в координатной форме, производится деление первого урав- нения системы на коэффициент при первом неизвестном, после че- 22 Перейти к оглавлению на странице: 26 го это преобразованное уравнение, будучи умноженным последо- вательно на все коэффициенты при первом неизвестном из всех нижележащих уравнениях системы, вычитается из соответствую- щих уравнений. Тем самым первое неизвестное оказывается “ис- ключенным” из всех уравнений системы, кроме первого. Оставив на время в покое первое уравнение, таким же образом исключа- ем из системы уравнений с номерами 2–n второе неизвестное и так далее. Результатом проделанной работы явится уравнение, содер- жащее лишь последнее неизвестное системы (2.1) , из которого оно и находится. Этот этап преобразований СЛАУ носит название прямо- го хода метода Гаусса. Вслед за найденным неизвестным находятся поочередно и остальные. Этот этап называется обратным ходом ме- тода Гаусса. Отметим очевидный недостаток простого метода Гаусса – воз- можное обращение в нуль ведущих элементов, т.е. тех элементов, на которые приходится делить в прямом ходе метода Гаусса. Легко указать и способ избавления от этого недостатка простого метода Гаусса – перестановка уравнений системы, которая, очевидно, не меняет решения СЛАУ. Ввиду того, что решается СЛАУ с неособой матрицей, при исключении первого неизвестного в первом столбце матрицы A найдётся ненулевой элемент и строка, содержащая этот элемент, переставляется на место первого уравнения. С целью ми- нимизации погрешности округления обычно выбирают не просто ненулевой элемент в столбце, а элемент, максимальный по модулю и этот процесс повторяется на каждом этапе прямого хода метода Гаусса. Данная модификация называется методом Гаусса с выбором максимального элемента по столбцу. Можно обобщить этот подход и выбирать максимальный элемент по всей матрице, что приводит к необходимости переобозначения неизвестных и потому использу- ется значительно реже. Отметим здесь, что если матрица A системы (2.1) обладает свойством диагонального преобладания (т.е. для всех строк мат- рицы A = {a ij } выполняются неравенства |a ii | > P j6=i |a ii | ), то обращение в нуль ведущих элементов в методе Гаусса не происхо- дит. В подробной записи метод Гаусса выглядит так: имеется 23 Перейти к оглавлению на странице: 26 СЛАУ, записанная в виде a 11 x 1 + a 12 x 2 + a 13 x 3 + . . . + a 1n x n = b 1 a 21 x 1 + a 22 x 2 + a 23 x 3 + . . . + a 2n x n = b 2 a 31 x 1 + a 32 x 2 + a 33 x 3 + . . . + a 3n x n = b 3 a n1 x 1 + a n2 x 2 + a n3 x 3 + . . . + a nn x n = b n (2.2) Пусть a 11 6= 0, иначе меняем порядок уравнений с тем, что- бы было выполнено данное условие. Разделим первое уравнение на a 11 , после чего оно примет вид: x 1 + a (1) 12 x 2 + a (1) 13 x 3 + . . . + a (1) 1n x n = b (1) 1 (2.3) где a (1) 1j = a 1j /a 11 , b (1) 1 = b 1 /a 11 Умножим уравнение (2.3) на a 21 и вычтем полученное уравне- ние из второго уравнения системы (2.2) . Аналогично преобразуем остальные уравнения. В итоге получим СЛАУ x 1 + a (1) 12 x 2 + a (1) 13 x 3 + . . . + a (1) 1n x n = b (1) 1 a (1) 22 x 2 + a (1) 23 x 3 + . . . + a (1) 2n x n = b (1) 2 a (1) n2 x 2 + a (1) n3 x 3 + . . . + a (1) nn x n = b (1) n , (2.4) где введены обозначения a (1) 2j = a 2j − a (1) 1j a 21 , a (1) 3j = a 3j − a (1) 1j a 31 , . . . , a (1) nj = a nj − a (1) 1j a n1 , b (1) j = b j − b (1) j a j1 , . . . , j = 2, n. Теперь, оставив без изменений первое уравнение системы (2.4) , можно применить описанную процедуру к системе из (n − 1) -го уравнений, исключив неизвестное x 2 из третьего и последующих уравнений. Получим систему вида x 1 + a (1) 12 x 2 + a (1) 13 x 3 + . . . + a (1) 1n x n = b (1) 1 x 2 + a (2) 23 x 3 + . . . + a (2) 2n x n = b (2) 2 a (2) n3 x 3 + . . . + a (2) nn x n = b (2) n , 24 Перейти к оглавлению на странице: 26 где a (2) 2j = a (1) 2j /a (1) 22 , b (2) 2 = b (1) 2 /a (1) 22 , a (2) 3j = a (1) 3j − a (2) 2j a (1) 32 , a (2) nj = a nj (1) − a (2) 2j a (1) n2 , b (2) j = b j (1) − b (2) j a (1) j , . . . , j = 3, n. Продолжая аналогичные вычисления, в итоге получим эквивалент- ную исходной треугольную систему x 1 + a (1) 12 x 2 + a (1) 13 x 3 + . . . + a (1) 1n x n = b (1) 1 x 2 + a (2) 23 x 3 + . . . + a (2) 2n x n = b (2) 2 x 3 + . . . + a (3) 3n x n = b (3) 3 x n = b (n) n , откуда и находятся последовательно все компоненты решения. ♣ 25 ОГЛАВЛЕНИЕ Курсовая работа по численным методам . . . . . . . . . . . . . . . 2 Глава 1. Численные методы решения задачи Коши для обыкновенных дифференциальных урав- нений . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 § 1. Постановка задачи . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 § 2. Одношаговые методы . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 Глава 2. Вычисление определенного интеграла . . . . . . . 10 § 1. Квадратурные формулы . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 § 2. Квадратурная формула Симпсона . . . . . . . . . . . . . . . . . . . . . . 12 Глава 3. Метод Ньютона решения уравнений . . . . . . . . . 17 § 1. Основы метода . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17 Глава 4. Задача аппроксимации . . . . . . . . . . . . . . . . . . . . . . . . 20 § 1. Линейная задача метода наименьших квадратов . . . . . . . . 20 § 2. Методы Гаусса . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 26 |