Д.П. Костомаров, А.П. Фаворский - Вводные лекции по численным методам (pdf) (1113729), страница 18
Текст из файла (страница 18)
Выше уже доказано, что для полиномов степени ( n − 1) она является точной.nПоследний переход заключается в том, что в сумму∑ c r ( x ) добавлены слагаемыеi =1i n −1iPn ( xi ) qn−1 ( xi ) . Они не меняют значения суммы, поскольку все равны нулю: ведьузлами квадратурной формулы являются корни полинома Лежандра Pn ( x ) .Итак, построенная квадратурная формула действительно является точной длялюбого полинома степени ( 2n − 1) , т.
е. задача Гаусса решена. На оценке погрешностиквадратурных формул Гаусса мы останавливаться не будем, однако задачи, к разборукоторых переходим, показывают, что эти формулы обеспечивают для гладкихфункций очень высокую точность.Задача 5.Построить квадратурную формулу Гаусса с двумя и тремя узлами.- 85 -Выведем сначала квадратурную формулу с двумя узлами. Узлы определяютсякак корни второго полинома Лежандра, выражение для которого мы выписываливыше (85). В данном случае имеем:(100)x1 = −1/ 3 , x2 = 1/ 3 .Узлы расположены симметрично относительно точки x = 0 .Весовые коэффициенты рассчитываются по формуле (95):11x − x2x − 1/ 3(101)c1 = ∫dx = ∫dx = 1 ,xx−−232−1 1−111x − x1x + 1/ 3c2 = ∫dx = ∫dx = 1 .xx−231−1 2−1Они равны между собой, а их сумма, в соответствии с общим соотношением (82),равна двум.
В результате искомая квадратурная формула принимает вид:1∫ f ( x ) dx = f ( −1/ 3 ) + f (1/ 3 ) + δ2.(102)−1Она является точной для любого полинома третьей степени.Перейдем теперь к выводу квадратурной формулы Гаусса с тремя узлами.Согласно формуле (85) для третьего полинома Лежандра ее узлами являются числа:(103)x1 = − 3 5 , x2 = 0 , x3 = 3 5 .Остается подсчитать весовые коэффициенты:(x x− 3 51c1 =∫−11c2 =∫)5dx = ,9− 3/ 5 −2 3/ 5()()( x + 3 5 )( x − 3 5 ) dx = 8 ,−3/ 5−1(x +1c3 = ∫−1)(104)935 x5dx = .93/ 5 2 3/ 5()()В результате квадратурная формула Гаусса с тремя узлами запишется в виде:15853/50fxdx=f−+f+f 3/ 5 + δ 3 .()()∫−1999Она является точной для любого полинома пятой степени.()()(105)Задача 6.Вычислить по формулам Симпсона и Гаусса при n = 2 интеграл:1∫ e dx = 2sh1 = 2.350402 .x−1Сравнить результаты численного интегрирования с точным значением интеграла имежду собой.- 86 -Формулы Симпсона и Гаусса дают в данном случае следующие результаты:14 2S 2 = ( e −1 + 4 + e ) = + ch1 = 2.362054 ,33 3γ 2 = -0.011651 ,1G2 = e −1/ 3 + e1/ 3 = 2ch= 2.342696 ,3δ 2 = 0.007706 .Мы видим, что даже с двумя узлами формула Гаусса дает хороший ответ.
Еготочность выше точности ответа, полученного по формуле Симпсона.В заключение сделаем следующее замечание. Несмотря на высокую точностьквадратурных формул Гаусса, при компьютерных расчетах ими пользуютсясравнительно редко. Дело в том, что для применения метода Гаусса нужно либоввести в компьютер до начала расчетов корни полинома Лежандра и весовыекоэффициенты, либо составить специальную подпрограмму для их вычисления. Врезультате потери человеческого и машинного времени на подготовку программы косновному расчету, связанному с вычислением интеграла, могут не окупитьсяточностью метода Гаусса.
Вычисление интеграла по более простой схеме методаСимпсона имеет с этой точки зрения преимущество.()§4. Построение первообразной с помощью численногоинтегрирования.Формулы Ньютона-Лейбница (1) позволяет выразить значение определенногоинтеграла от функции f ( x) через ее первообразную F ( x) .
В математическом анализеустанавливается и прямо противоположная возможность: первообразная функцииf ( x) , непрерывной на отрезке [ a, b ] , может быть записана в виде определенногоинтеграла с переменным верхним пределом:xF ( x ) = ∫ f ( t ) dt .(106)x0Здесь x0 , x - две точки отрезка [ a, b ] , причем нижний предел интегрирования x0предполагается фиксированным, верхний x - переменным. В случае непрерывнойфункции f ( x) функция F ( x) , определенная с помощью интеграла (106), являетсядифференцируемой и ее производная равна f ( x) :x⎞d ⎛F ′ ( x ) = ⎜ ∫ f ( t ) dt ⎟ = f ( x ) .⎟dx ⎜⎝ x0⎠(107)Формула (106) в сочетании с какой-нибудь формулой численного интегрирования,например, Симпсона, представляет собой универсальный алгоритм построенияпервообразной.
Приведем два примера, иллюстрирующие этот алгоритм.Функция f ( x ) = sin x/x непрерывна и, следовательно, имеет первообразные. Онине могут быть выражены через элементарные функции, но представление в видеинтеграла с переменным верхним пределом для них справедливо. Одну из- 87 -первообразных мы получим, выбирая нижний предел интегрирования x0 = 0 . Ееназывают интегральным синусом и обозначаютxsin tSi ( x ) = ∫dt .t0Интегральный синус определен на всей числовой прямой, является нечетнойфункцией x , имеет конечные предельные значения на бесконечностиlim Si ( x ) = ±x →±∞Согласно (107)π2.Si′ ( x ) = sin x / x .По знаку производной легко определить области возрастания и убывания функции,разделенные точками экстремума xk = kπ ( k = ±1, ± 2, K) .
Методы численногоинтегрирования позволяют вычислить значения Si ( x ) при любом x . Графикинтегрального синуса при x ≥ 0 приведен на рис. 3.В качестве второго примера рассмотрим функцию ошибок erf ( x ) , играющуюважную роль в теории вероятности. Ее обозначение образовано с помощью первыхбукв английского названия функции ошибок – error function. Подобно интегральномусинусу, функция ошибок вводится в виде интеграла с переменным пределом от2функции e − x , которая не имеет первообразных в классе элементарных функций:x2−t 2erf ( x ) =e∫ dt .π 0Функция ошибок определена на всей числовой прямой, является нечетной функциейx , имеет конечные предельные значения на бесконечности:lim erf ( x ) = ±1 .x →±∞Согласно (107)erf ′ ( x ) =22e− x .πПроизводная всюду положительная, следовательно, функция ошибок монотонновозрастает. Ее график приведен на рис.
4.Существует ряд других специальных функций, которые вводятся как интегралыс переменным верхним пределом. Не будем останавливаться на их описании, отметимлишь, что разобранные примеры показывают, насколько условно деление функций наэлементарные и неэлементарные. По существу, чтобы работать с какой-нибудьфункцией, нужно знать ее свойства и иметь алгоритм вычисления при любом значенииаргумента. С этой точки зрения применение интегрального синуса или функцииошибок ничем не отличается от применения привычных нам элементарных функций.- 87 -Глава 5. ЧИСЛЕННОЕ ИНТЕГРИРОВАНИЕ ОБЫКНОВЕНЫХДИФФЕРЕНЦИАЛЬНЫХ УРАВНЕНИЙНаиболее универсальными методами численного решения обыкновенныхдифференциальных уравнений являются разностные методы. Они основаны на заменепроизводных в дифференциальном уравнении разностными отношениями.
Врезультате исходное дифференциальное уравнение сводится к системе алгебраическихуравнений, которые называются разностными. Решение этой системы даетприближенное решение исходной задачи.§1. Разностная аппроксимация производных.1.1. Сеточные функции.Пусть на отрезке [ a, b ] задан набор точекa = x0 < x1 < x2 < L < xn−1 < xn < b .(1)Будем называть его сеткой. Чтобы не усложнять изложения, условимся считать сеткуравномерной:(2)xi +1 − xi = h = ( b − a ) / n , 0 ≤ i ≤ n − 1 .Пусть каждой точке сетки xi сопоставлено по определенному закону число yi .Совокупность этих чисел y = { y0 , y1 ,K , yn } = { yi } ( 0 ≤ i ≤ n ) назовем сеточнойфункцией. Сеточные функции, определенные на сетке (1), образуют ( n + 1) -мерноелинейное пространство.Чтобы иметь возможность сравнивать сеточные функции между собой, говорить обих близости, нужно ввести в этом пространстве норму.
В этой главе мы будемпользоваться нормой С , которая определяется следующим образом:y c = max yi .(3)0≤i ≤nЭто определение законно, поскольку удовлетворяет трем аксиомам нормы:1. Норма неотрицательнаy c ≥ 0,причем равенство нулю имеет место только для нулевого элемента.2. Модуль числового множителя можно вынести за знак нормыαy c = α y .3. Неравенство треугольникаy+z c ≤ y c + z c.Справедливость последнего утверждения вытекает из свойства максимума:max yi + zi ≤ max yi + max zi .0≤i ≤ n0≤ i ≤ n(4)0 ≤i ≤ n1.2.
Разностные аппроксимации первой производной.Для сеточных функций нельзя ввести обычное понятие производной,включающее операцию предельного перехода при ∆x → 0 . Вместо производной здесьвводятся разностные отношения:- 88 -yi +1 − yi, 0 ≤ i ≤ n −1;(5)hy − yi −1L−h [ yi ] = i,1 ≤ i ≤ n ;(6)hyi +1 − yi −1L(0),1 ≤ i ≤ n − 1 .(7)h [ yi ] =2hОтношение (5) называют правой разностной производной, отношение (6) – левойразностной производной и отношение (7) – центральной разностной производной.Чтобы установить связь разностных отношений (5) – (7) с обычной производной,предположим, что на отрезке [ a, b ] определена дифференцируемая функция y ( x ) ,значения которой в точках сетки (1) равны значениям рассматриваемой сеточнойфункции: yi = y ( xi ) .
Вычислим первую производную функции y ( x ) в точках xi исопоставим с разностными отношениями (5) – (7):(8)ψ i+ = L+h [ yi ] − y′ ( xi ) , 0 ≤ i ≤ n − 1 ;L+h [ yi ] =(9)ψ i− = L−h [ yi ] − y′ ( xi ) ,1 ≤ i ≤ n ;′(10)ψ i(0) = L(0)h [ yi ] − y ( xi ) ,1 ≤ i ≤ n − 1 .Эти величины представляют собой погрешности аппроксимации производной спомощью разностных отношений (5) – (7) в точке xi .Предположим, что функция y ( x ) дважды непрерывно дифференцируема наотрезке [ a, b ] и запишем для нее формулу Тейлора с остаточным членом в формеЛагранжа1yi +1 = y ( xi + h ) = yi + y ′ ( xi ) h + y ′′ ( xi + θ i h ) h 2 ,(11)2где θ i какое-то неизвестное нам число между нулем и единицей. Подставляяразложение (11) в формулу (8), получим1ψ i+ = y ′′ ( xi + θi h ) h .(12)2Аналогичное представление можно получить для величины ψ i− (9)1ψ i− = − y ′′ ( xi − θ i h ) h .(13)2Формулы (12) и (13) не позволяют вычислить соответствующие погрешности, но даютвозможность их оценить.
Функция y′′ ( x ) , по предположению, непрерывна на отрезке[ a, b] , и, следовательно, ограничена:(14)y′′ ( x ) ≤ M 2 , a ≤ x ≤ b .В результате получаем11ψ i+ ≤ M 2 h , ψ i− ≤ M 2 h .(15)22- 89 -Оценки (15) являются равномерными, поскольку не зависят от индекса i . Такимобразом, левое и правое разностное отношение аппроксимируют производную y′ ( x ) спервым порядком точности относительно h .Для оценки ψ i(0) (10) предположим, что функция y ( x ) три раза непрерывнодифференцируема на отрезке [ a, b ] и продолжим разложение (11) еще на один член1yi +1 = yi + y ′ ( xi ) h + 12 y ′′ ( xi ) h 2 + y ′′′ ( xi + θ1,i h ) h 3 ,6(16)1231yi −1 = yi − y ′ ( xi ) h + 2 y ′′ ( xi ) h − y ′′′ ( xi − θ 2,i h ) h .6Подставляя разложения (16) в формулу (10), будем иметь1(17)ψ i(0) = y ′′′ ( xi + θ1,i h ) + y ′′′ ( xi − θ 2,i h ) h 2 .6По предположению функция y′′′ ( x ) непрерывна и, следовательно, ограничена наотрезке [ a, b ] :(18)y′′′ ( x ) ≤ M 3 , a ≤ x ≤ b .{}В результате из равенства (17) получим оценку1ψ i(0) ≤ M 3h 2 .(19)3Оценка (19), как и раньше (15), не зависит от индекса i , она является равномерной.Таким образом, центральная разностная производная дает более хороший результат:она аппроксимирует производную y′ ( xi ) со вторым порядком точности относительноh для функций, трижды непрерывно дифференцируемых на отрезке [ a, b ] .Задача 1.Рассмотреть функцию y = 1/ (1 − x ) на сеткеx0 = −0.1, x1 = 0 , x2 = 0.1 .(20)Вычислить в точке x1 = 0 правую, левую и центральную разностные производные,найти погрешности аппроксимации производной y′(0) = 1 , сравнить их с априорнымиоценками по формулам (15) и (19).В данном случае1/ 0.9 − 1= 1.111111,0.1ψ 1+ = 0.111111 ;1 − 1/1.1L−h [ y1 ] == 0.909090 ,0.1ψ i− = −0.090909 ;1/ 0.9 − 1/1.1L(0)= 1.010101 ,h [ y1 ] =0.2L+h [ y1 ] =- 90 -ψ i(0) = 0.010101 .Перейдем к априорной оценке погрешности.