Интерполяция функций одной и нескольких переменных

Содержание Лекция_2 ИНТЕРПОЛЯЦИЯ ФУНКЦИЙ ОДНОЙ И НЕСКОЛЬКИХ ПЕРЕМЕННЫХ Лекция_3 Нелинейная интерполяция Лекция_4 НАИЛУЧШЕЕ СРЕДНЕКВАДРАТИЧНОЕ ПРИБЛИЖЕНИЕ Лекция_5 ЧИСЛЕННОЕ ИНТЕГРИРОВАНИЕ ФУНКЦИЙ ОДНОЙ ПЕРЕМЕННОЙ Лекция_6 КРАТНЫЕ ИНТЕГРАЛЫ Лекция_7 ЧИСЛЕННОЕ ДИФФЕРЕНЦИРОВАНИЕ Лекция_8 ЧИСЛЕННОЕ ДИФФЕРЕНЦИРОВАНИЕ Лекция_9-10 ОБЫКНОВЕННЫЕ ДИФФЕРЕНЦИАЛЬНЫЕ УРАВНЕНИЯ Лекция_11 МОДЕЛИ НА ОСНОВЕ ОДУ. КРАЕВАЯ ЗАДАЧА Лекция_12 МОДЕЛИ НА ОСНОВЕ ОДУ. КРАЕВАЯ ЗАДАЧА Лекция_13 МОДЕЛИ НА ОСНОВЕ ОДУ. КРАЕВАЯ ЗАДАЧА Лекция_14-17 МОДЕЛИ НА ОСНОВЕ ДИФФЕРЕНЦИАЛЬНЫХ УРАВНЕНИЙ В ЧАСТНЫХ ПРОИЗВОДНЫХ (ДУЧП) 2 8 13 17 26 31 38 42 55 60 66 71 1. ИНТЕРПОЛЯЦИЯ ФУНКЦИЙ ОДНОЙ И НЕСКОЛЬКИХ ПЕРЕМЕННЫХ В задачах математического моделирования, которые для своей реализации требуют, как правило, применения ЭВМ, из различных способов задания функций наиболее часто используется табличный способ. При этом возникает необходимость определения значений функции для значений аргумента, отличающихся от тех, которые содержатся в таблице. В этих целях применяют процедуру интерполяции, а также экстраполяции. По сути, речь идет о восстановлении непрерывной функции по заданному набору ее значений для счетного множества значений аргумента. Задача восстановления функции может ставиться и иначе, когда ищется функция, близкая в определенном смысле к таблично заданной функции (например, наилучшее среднеквадратичное приближение). В обоих случаях исходную функциональную зависимость y( x ) , по которой сопоставлена таблица, заменяют приближенной функцией  ( x , a ) , которую и используют в дальнейших вычисления. При этом близость y( x ) и  ( x , a ) обеспечивается подбором свободных параметров a  a0 , a1 ,..., a n . В простейшем варианте лагранжевой интерполяции, когда функция y( x ) задана в узлах некоторой сетки x i , параметры a определяют из условия совпадения  ( x , a ) со зна- чениями функции в фиксированном числе узлов. Получаемая при этом система уравнений имеет вид ( xi , a0 , a1 ,..., a n )  y( xi )  yi ,0  i  n. (2.1) Из этой системы можно определить все ai . Интерполяция может быть линейной или нелинейной, в соответствии с характером зависимости функции  ( x , a ) от параметров a . 1.1. Линейная интерполяция При линейной интерполяции функция  ( x , a ) имеет вид n ( x , a )   a j  j ( x ), j 0 (2.2) где все функции  j ( x ) линейно независимы. Подставляя (2.2) в (2.1), получим систему линейных уравнений для определения ai : n a  j 0 j j ( xi )  y i , 0  i  n. (2.3) Единственность решения обеспечивается требованием неравенства нулю определителя системы (1.3):  0 ( x0 )  1 ( x0 )...  n ( x0 )  . ... . 0 (2.4)  0 ( x n )  1 ( x n )...  n ( x n ) при x k  xl (т.е. при любых несовпадающих узлах). Систему функций, удовлетворяющую условию (2.4), называют чебышевской. Из различных систем функций наиболее распространены многочлены, хотя применяют также тригонометрические и экспоненциальные функции. Если в качестве системы функций выбрать степенные функции аргумента, т.е.  j ( x )  x j , то определитель (2.4) окажется определителем Вандермонда, который не равен нулю при условии x k  xl . Следовательно, интерполяционный полином всегда существует и он единствен. Для практических вычислений удобно использовать многочлен Pn ( x ) в форме интерполяционного полинома Ньютона. Введем понятие разделенных разностей функции y( x ) , заданной в узлах x i : y( x i , x j )  [ y( x i )  y( x j )] /( x i  x j ), y( x i , x j , x k )  [ y( x i , x j )  y( x j , x k )] /( x i  x k ), (2.5) y( x i , x j , x k , x l )  [ y( x i , x j , x k )  y( x j , x k , x l )] /( x i  x l ), Правило образования таких конструкций понятно из приведенной записи. Разделенные разности первого, второго и более высоких порядков связаны с производными соответствующих порядков. 2 Полином n-ой степени может быть представлен через разделенные разности в виде: Pn ( x )  Pn ( x0 )  ( x  x0 )P( x0 , x1 )  ( x  x0 )( x  x1 )P( x0 , x1 , x 2 ) ...  ( x  x0 )( x  x1 )...( x  x n 1 )P( x0 , x1 , x 2 ,..., x n ) (2.6) Ввиду того, что значения интерполяционного полинома в этих узлах совпадают со значениями интерполируемой функции, разделенные разности, входящие в (2.6), выражаются через узловые начения функции. В результате получается полином, называемый интерполяционным многочленом Ньютона: n y( x )  y( x0 )   ( x  x0 )( x  x1 )...( x  x k 1 ) y( x0 ,..., x k ). (2.7) k 1 При вычислениях по этой формуле точность расчетов удобно оценивать, наблюдая за тем, насколько быстро убывают члены ряда. Если это происходит достаточно быстро, можно оставлять только те члены, которые больше заданной погрешности расчетов. Погрешность многочлена Ньютона можно оценить по формуле y( x )  Pn ( x )  где M n 1  max y  n 1 M n 1  n( x ), ( n  1 )! (  ) - максимальное значение производной интерполируемой функ- ции на отрезке между наименьшим и наибольшим из значений x0 , x1 , x 2 ,..., x n , а полином n  n ( x )   ( x  xi ) . i 0 На равномерной сетке узлов приведенная формула трансформируется к виду y( x )  Pn ( x )  h M n 1    2 2n n 1 , (2.8) ,где h  xi 1  xi  xi  xi 1     (предполагается, что шаг сетки постоянен). Выражение (2.8) свидетельствует, что с уменьшением расстояния между узлами (шага) h погрешность представления функций полиномом Ньютона убывает как O( h n 1 ). Трудность использования указанных теоретических оценок на практике состоит в том, что производные интерполируемой функции обычно неизвестны, поэтому для определения погрешности удобнее воспользоваться оценкой первого отброшенного члена. Пример. Построить интерполяционный полином Ньютона четвертой степени для функции y( x )  cos( 2 x ) в области значений аргумента 0  x  / 4 . Заполним таблицу разделенных разностей, вычисляемых по пяти узлам, представив для удобства вычислений y( z )  cos[(  / 2 ) z ] , 0  z  1. z0 y( z 0 ) y( z 0 , z 1 ) y( z 0 , z 1 , z 2 ) y( z 0 ,...,z 3 ) y( z 0 ,...,z 4 ) 1 -0,304 0,25 0,924 -1,128 -0,868 0,5 0,707 0,363 -0,856 -1,296 0,75 0,383 0,149 0,512 -0,472 -1,532 1 Полином Ньютона y ( z )  1  0,304 z  1,128 z ( z  0,25)  0,363z ( z  0,25)( z  0,5)   0,149 z ( z  0,25)( z  0,5)( z  0,75) Вычислим y( 0 ,6 )  1  0 ,182  0 ,237  7 ,623  10 3  4 ,694  10 3  0 ,589. Точное значение y( 0,6 )  0,588 , т.е. точность вычислений по приближенной формуле оказалась весьма высокой. Опишем алгоритм построения полинома Ньютона. В качестве исходных данных задаются: таблично интерполируемая функция, значение аргумента x, для которого выполняется интерполяция, и степень полинома. 1. Формируется конфигурация из (n+1)-го узла, по возможности симметрично относительно значения x (понятно, что вблизи краев таблицы это невозможно). 2. На сформированной конфигурации узлов строится таблица, аналогичная той, которая представлены в вышеприведенном примере. 3. Строится полином Ньютона с разделенными разностями из верхней строки таблицы. Обратная интерполяция. Применяется для приближенного нахождения корня монотонных функций. Суть ее состоит в том, что меняются местами столбцы табличной функции и выполняется обычная интерполяция при задании x=0. Факт наличия корня у функции устанавливается по наличию смены знака у функции при продвижении по строкам таблицы. Помимо полинома Ньютона существует еще одна важная форма интерполяционного многочлена- полином в форме Лагранжа n Ln ( x )   L(i n ) ( x ) y( xi ) , i 0 (n) где лагранжевы коэффициенты Li ( x ) определяются по формуле ( x  xk ) . k 0 ( x i  x k ) n L(i n ) ( x )   k i Полином Лагранжа в отличие от полинома Ньютона содержит в явном виде значения интерполируемой функции в узлах y( x i ) . В заключение отметим, что в практике вычислений для интерполяции полиномы степени выше пятой обычно не используют, т.е. число узлов интерполяции не превышает шести. Если при таком числе узлов не обеспечивается заданная погрешность, следует уменьшать расстояние между узлами. 1.2. Многомерная интерполяция В различных приложениях широко используют двумерные и трехмерные таблицы. Например, теплофизические свойства различных веществ зависят от температуры и давления, а оптические характеристики - еще и от длины волны излучения. Осуществляя многомерную интерполяцию, следует помнить, что расположение узлов не может быть произвольным. Например, при интерполяции полиномом первой степени P1 ( x , y ) узлы не должны лежать на одной прямой в плоскости.. Проверять условия подобного типа достаточно сложно, поэтому на практике целесообразно строить регулярные сетки, как правило, прямоугольные и равномерные, когда узлы являются точками пересечения двух взаимно перпендикулярных систем параллельных прямых. На этой сетке проводят простую последовательную интерполяцию: сначала по строкам, а затем по столбцам. 1.3. Нелинейная интерполяция Для табулирования быстроменяющихся функций требуется весьма малый шаг, т.е. возникает необходимость создавать таблицы очень больших объемов, что в ряде случаев неприемлемо. Оказывается, что преобразованием переменных   ( y ) и добиться того, чтобы в новых переменных график    ( x ) можно  (  ) был близок к прямой хотя бы на отдельных участках. В этом случае интерполяцию проводят в переменных ( , ) , а затем обратным интерполированием находят y i  y(  i ) . Преобразования  ( y ) и  ( x ) должны быть достаточно простыми (логарифмиче- ская, экспоненциальная, тригонометрические и некоторые другие функции). При этом надо заботиться о том, чтобы и обратное преобразование y(  ) оказалось несложным. Во многих задачах теплофизики, гидродинамики, оптики, переноса излучения) и других областей науки и техники часто встречается степенная зависимость функции от своих аргументов. В этом случае удобны преобразования типа логарифмирования. Вообще говоря, в каждом конкретном случае приходится специально подбирать вид функций  ( y )и   ( x ) . Например. функция y  a может быть представb  c  e dx лена следующим образом ln( a a  b )  d  x  ln c , т.е.   ln(  b ),   x . y y Для функциональной зависимости y  переменных   x линеаризация достигается заменой ax  b y 1 1 ,   , а для функции y  axe  bx -   ln( ),   x . x y x Пример. Получить формулу для нелинейной двухточечной интерполяции функции   ln( y ) и y( x ) , если переменные можно преобразовать по формулам   x. Составим интерполяционный полином Ньютона на двухточечном шаблоне:   0  1  0 (    0 ). 1  0 В исходных переменных имеем ln( y )  ln( y 0 )  ln( y 1 )  ln( y 0 ) ( x  x 0 ), x1  x0 и окончательно y  y 0 ( y 1 / y 0 )( x  x ) /( x1  x0 ) . 1.4. Интерполяция сплайнами Слово “сплайн” переводится как “гибкая линейка”. Такую линейку можно использовать для проведения кривых через заданную совокупность точек, изгибая и придерживая ее так, чтобы ребро проходило через все точки на плоскости. Равновесие гибкой линейки описывается уравнением  IV ( x )  0 , т.е. интерполяционный полином на участке между каж- дой парой соседних точек имеет третью степень:  ( x )  a i  bi ( x  xi 1 )  ci ( x  xi 1 ) 2  d i ( x  xi 1 ) 3 , (3.1) xi 1  x  xi ,0  i  N . В узлах значения многочлена и интерполируемой функции совпадают:  ( x i 1 )  y i 1  a i , (3.2)  ( xi )  y i  ai  bi hi  ci hi2  d i hi3 , (3.3) hi  xi  xi 1 ,1  i  N . Число таких уравнений меньше числа неизвестных в два раза. Недостающие уравнения получают, приравнивая во внутренних узлах первые и вторые производные, вычисляемые по коэффициентам на соседних участках:  ( x )  bi  2ci ( x  xi 1 )  3d i ( x  xi 1 ) 2 ,  ( x )  2ci  6 d i ( x  xi 1 ), xi 1  x  xi , bi 1  bi  2ci hi  3d i hi2 , (3.4) ci 1  ci  3d i hi , 1  i  N  1. (3.5) Недостающие условия можно получить, полагая, например, что вторая производная равна нулю на концах участка интерполирования:  ( x0 )  0,c1  0, (3.6)  ( x N )  0,c N  3d N hN  0, (3.7) Уравнения (3.2)-(3.7) позволяют определить все 4 N неизвестных коэффициентов: ai ,bi , ci , d i ( 1  i  N ). Решение полученной системы уравнений можно сильно упростить, если привести ее к специальному виду. Используя уравнение (3.2), можно получить все коэффициенты ai . Из (3.5) и (3.7) следует di  c i 1  c i , 3hi dN   1  i  N  1. (3.8) cN 3h N (3.9) Из (3.3) и (3.8): bi  y i  y i 1 c  2c i  hi i 1 , hi 3 Из (3.3) и (3.9): 1 i  N 1 . (3.10) bN  y N  y N 1 2c  hN N hN 3 (3.11) Исключим теперь из (3.4) величины bi и bi 1 с учетом (3.10), наращивая во втором случае индекс на 1, а величину d i - с учетом (3.8). В результате получим систему уравнений для определения коэффициентов c i : c1  0 , hi 1 c i 1  2( hi 1  hi )c i  hi c i 1  3( y i  y i 1 y i 1  y i  2  ) , 2  i  N  1 (3.12) hi hi 1 c N 1  0 . ) После нахождения коэффициентов c i остальные коэффициенты определяют по следующим формулам: a i  y i 1 , 1  i  N , bi  y i  y i 1 c  2c i  hi i 1 , 1 i  N 1 , hi 3 bN  y N  y N 1 2c  hN N , hN 3 di  ci 1  ci , 1 i  N 1 , 3hi dN   cN . 3h N Осталось выяснить, как решать систему (3.12). Матрица этой системы трехдиагональна, т.е. все ее элементы равны нулю, кроме тех, которые находятся на главной и двух соседних диагоналях. Такие системы удобно решать методом прогонки. Суть метода в следующем. Применение метода исключения Гаусса для решения системы уравнений с трехдиагональной матрицей приводит к тому, что система уравнений преобразуется к виду, когда в каждом уравнении содержится только два неизвестных, и при обратном ходе одно из этих неизвестных выражается через другое. Поэтому применительно к (3.12) можно записать: ci   i  1 ci  1   i  1 , где (3.13)  i 1 , i 1 - некоторые не известные пока прогоночные коэффициенты; c i 1   i c i   i . Подставляя последнее выражение в (3.12) и преобразуя, получим ci   hi f i  hi 1 i ci 1  . hi 1 i  2( hi 1  hi ) hi 1 i  2( hi 1  hi ) (3.14) Сравнивая (3.13) и (3.14), получим  i 1   hi f i  hi 1 i , i 1  . hi 1 i  2( hi 1  hi ) hi 1 i  2( hi 1  hi ) (3.15) В этих формулах введено обозначение f i  3( y i  y i 1 y i 1  y i  2  ). hi hi 1 Из условия c1  0 , следует  2  0 , 2  0 . Теперь алгоритм решения (3.12) выглядит следующим образом. По формулам (3.15) при известных  2 , 2 , равных нулю, вычисляют прогоночные коэффициенты  i 1 , i 1 ( 2  i  N ) (прямой ход). Затем по формулам (3.13) при условии c N 1  0 определяют все c i (обратный ход). ЛЕКЦИЯ №4. НАИЛУЧШЕЕ СРЕДНЕКВАДРАТИЧНОЕ ПРИБЛИЖЕНИЕ Процедура интерполяции функции подразумевает построение некоторой новой функции, совпадающей с заданной в фиксированных узлах. В ряде случаев целесообразно приближать функции не по точкам, а в среднем, например, когда значения функции в узлах определены неточно. Пусть имеется множество функций  ( x ) , принадлежащих линейному пространству функций. Под близостью в среднем исходной y и аппроксимирующей  функций будем понимать результат оценки суммы N I   i [ y ( xi )   ( xi )]2 (4.1) i 1 где  i - вес точки. Суммирование выполняется по всем N узлам заданной функции. Такой вид аппроксимации называют среднеквадратичным приближением. Можно рассмотреть две задачи: 1 - подобрать функцию  ( x ) так, чтобы выполнялось неравенство I  ; 2 - найти наилучшее приближение, т.е. такую функцию  ( x ) , чтобы было справед- ливым соотношение N  i 1 i [ y( x i )  ( x i )] 2  min (4.2) Далее займемся отысканием наилучшего приближения, которое применительно к таблично заданным функциям называется методом наименьших квадратов. Разложим функцию  ( x ) по системе линейно независимых функций  k ( x ) : n  ( x )   ak  k ( x ) . (4.3) k 0 В дальнейшем для сокращения записи будем пользоваться определением скалярного произведения в пространстве дискретно заданных функций N ( f , )   i f ( xi )  ( xi ), i  0 . i 1 Несложно установить, что имеют место следующие равенства, справедливые для обычного скалярного произведения элементов линейного пространства 1. ( f , )  ( , f ) 2. ( f   ,  )  ( f ,  )  ( ,  ) (4.4) Подставляя (4.3) в условие (4.2), получим с учетом (4.4.) n n n (( y   ), ( y   ))  ( y, y )  2 ak ( y, k )   ak am ( k , m )  min . k 0 k 0 m0 Дифференцируя это выражение по a k и приравнивая производные нулю, найдем n  ( , m0 k m ) a m  ( y ,  k ), 0k n . (4.5) Определитель этой системы в силу линейной независимости функций  k (x ) не равен нулю. Следовательно, из системы (4.5) можно найти коэффициенты a k , определяющие функцию  ( x ) согласно (4.3) и минимизирующие (4.1). Таким образом, наилучшее средне- квадратичное приближение существует и оно единственно. В качестве  k ( x ) чаще всего используют полиномы Лежандра, Чебышева, Лагерра, Эрмита, ортогональные с заданным весом. Наиболее употребительный вариант метода наименьших квадратов соответствует случаю степенного вида функций  k ( x ) , т.е.  k ( x )  x k , причем 0  k  n . Обычно в сумме (4.3) берут не более пяти-шести членов. Система уравнений (4.5) при этом принимает вид n  (x k m 0 где , x m ) a m  ( y, x k ) , 0  k  n , N (4.6) N ( x k , x m )    i xik  m , ( y, x k )    i y i xik . i 1 i 1 Пример 2.1. Методом наименьших квадратов аппроксимировать функцию линейной зави- ( x )  a 0  a 1 x . симостью вида В данном случае n  1. В итоге система уравнений (4.6) имеет вид ( x 0 , x 0 ) a0  ( x 0 , x 1 ) a1  ( y , x 0 ), ( x 1 , x 0 ) a0  ( x 1 , x 1 ) a1  ( y , x 1 ). Скалярные произведения в полученной системе записываются следующим образом: N N N ( x , x )    i , ( x , x )    i x i , ( x , x )    i xi2 1 1 i 1 1 i 1 i 1 N N i 1 i 1 ( y , x 0 )    i y i , ( y , x 1 )    i y i xi . Окончательно N a0  N i 1 i 1 N N i 1 N   x i i 1 N a1  N   i yi   i xi2    i xi   i xi yi i i 1 N   x y i i 1 i i 1 N i i i 1 i ,  (   i xi ) 2 i 1 N    i xi   i y i i 1 N i 2 i N   x i 1 i 1 N i 1 . N 2 i  (   i xi ) 2 i 1 k Система функций x не ортогональна, поэтому при больших n задача (4.5) плохо обусловлена, в связи с чем на практике ограничиваются значениями n  5 . На рис.1 представлен результат применения метода наименьших квадратов для аппроксимации данных (черные точки) полиномами степеней 1, 2, 4 и 6. Веса всех точек приняты равными 1. Рис. 1. Аппроксимация функции методом наименьших квадратов. Несколько замечаний о введенных выше (4.1) весах точек  i . Чем больше вес точки, тем ближе к точке проходит аппроксимирующая кривая. Под весом можно понимать, например, величину, обратную относительной погрешности задания функции, т.е. чем более точное значение имеет табличная функция в некоторой точке, тем больше ее вес и тем ближе к ней пройдет график аппроксимирующей функции. Подведем итоги. Для применения метода наименьших квадратов в случае аппроксимации полиномом следует действовать следующим образом. 1. Выбирается степень полинома n<0 решение существует, единственно, и прогонка устойчива в силу преобладания диагонального элемента матрицы системы: из (15), (16) ясно, что модуль этого элемента больше суммы модулей недиагональных членов. Видно, что при  =0 схема (9) переходит в явную схему (8), а при  =1- в чисто неявную (7). При   1 схема (9) называется симметричной (по времени). 2 ЛЕКЦИЯ №16. МОДЕЛИ НА ОСНОВЕ ДУЧП. ОСНОВНЫЕ ПОНЯТИЯ ОБ АППРОКСИМАЦИИ РАЗНОСТНЫХ СХЕМ Качество построенной разностной схемы оценивается такими свойствами как аппроксимация, устойчивость, сходимость. Ниже будет показано, что из аппроксимации и устойчивости разностной схемы следует сходимость приближенного решения к точному. Рассмотрим последовательно указанные свойства схем. Введем понятие невязки разностной схемы, построенной для дифференциального уравнения, записанного в общем операторном виде Au( x)  f , (17) с дополнительными условиями Bu ( x)   ( x) , (18) где, например, для уравнения (1) оператор A имеет вид  2  A    a 2  . x   t Разностная схема для задачи (17), (18): Ah y  h , (19) Bh y   h (20) Если подставить в соотношения (19) точное решение, то данное равенство будет нарушено, так как приближенное решение y не совпадает с точным решением u . Невязкой называется величина   h  Ah u  ( Au  f ) ( Ah u  h ) . (21) Для граничных условий получаем невязку в виде    h  Bh u  ( Bu   ) ( Bh u   h ) . (22) Дадим определение аппроксимации. Разностная схема (19), (20) аппроксимирует задачу (17), (18), если в некоторой норме   0,   0 при h  0 , и аппроксимация имеет   O( h p ) , p-ый порядок, если   O(h p ) при h  0 . Фигурирующие в приведенных соотношениях нормы могут быть определены как сеточные аналоги различных норм: чебышевской u ( x) C , гильбертовой u ( x) L , 2 энергетической u ( x ) E : u( x ) C  max  u( x ) , a  x b b   ( x) u u( x ) L  2 2 ( x) dx ,  ( x)  0 , a b  [  ( x) u u( x ) E  1 2 x ( x)   0 ( x) u 2 ( x)] dx , 1 ( x)  0,  0 ( x)  0 . a Указанные сеточные аналоги выписанных норм представляют в таком виде, чтобы при h  0 они переходили в эти нормы: y C  max  yn , 0 n  N y L2  N  n 1 n y n2 hn . Невязку оценивают, проводя разложение точного решения в ряд Тейлора. Найдем невязку разностной схемы (9) для уравнения (1). Выполним разложение решения на сетке, принимая за центр разложения точку  ( x n , t m  ) . Получим 2  1   1   1   h2 h3 h4 un 1  u  u t    u tt    u ttt    u tttt  hu x  u xx  u xxx  u xxxx  2 22 62 24  2  2 6 24 2 h 3 4 1   1   1   1   1    u tx    huttx   h 2 u txx    hutttx    h 2 u ttxx   h 3u txxx  ... 2 22 22 62 42 62 2 3 2 ,  1   1   1   un  u  u t    u tt    u ttt    u tttt  ... , 2 2!  2  3!  2  4!  2  2 u n 1  3 4 1   1   1   h2 h3 h4  u  u t    u tt    u ttt    u tttt  hu x  u xx  u xxx  u xxxx  2 22 62 24  2  2 6 24 h 2 3 4 1   1   1   1   1    u tx    huttx   h 2 u txx    hutttx    h 2 u ttxx   h 3u txxx  ... 2 22 22 62 42 62 u n 1 2 3  2 1   1   1   h2 h3 h4  u  u t    u tt    u ttt    u tttt  hu x  u xx  u xxx  u xxxx  2 22 62 24  2  2 6 24 h 2 3 4 1   1   1   1   1    u tx    huttx   h 2 u txx    hutttx    h 2 u ttxx   h 3u txxx  ... 2 22 22 62 42 62 2 3 2  1   1   1   u n  u  u t    u tt    u ttt    u tttt  ... 2 2!  2  3!  2  4!  2  2 3 4 В этих формулах введены обозначения: ut  u  2u u , u tt  2 , u x  и т.д. t x t , , Подставив эти разложения в формулу для невязки, придем к соотношению t t m   u ( x, t )   2 u ( x, t )    a  f ( x, t )  2 x  t  x  xn  2   un  un   a    u n 1  2u n  u n 1  h2 u n 1  2u n  u n 1  n  h2 1 2 1 ah 2   a (  ) u txx  (au ttxx  u ttt )  u xxxx   n  f ( x n , t m  )  O( 2  h 2 ). 2 8 3 12 2  a(1   ) Можно заметить, что если взять   n  f ( xn , t m  ) , 2 то при рассматриваемая разностная схема имеет аппроксимацию O( 2  h 2 ) , а при     1 2 1 2 O(  h 2 ) . Аналогично проверяется аппроксимация начальных и граничных условий, если они содержат производные от функции, например, граничных условий второго или третьего рода. ЛЕКЦИЯ №17. МОДЕЛИ НА ОСНОВЕ ДУЧП. ОСНОВНЫЕ ПОНЯТИЯ ОБ УСТОЙЧИВОСТИ И СХОДИМОСТИ РАЗНОСТНЫХ СХЕМ 17.1. Устойчивость Под устойчивостью задачи понимают непрерывную зависимость решения от входных данных, т. е. малые отклонения во входных данных должны приводить к малому изменению решения. Неустойчивость проявляется в том, что малые ошибки, допущенные на любом расчетном шаге, приводят к быстрому их нарастанию в ходе дальнейших вычислений, что, естественно, обесценивает получаемые результаты. Дадим определение устойчивости. Разностная схема Ah y  h , Bh y   h устойчива, если ее решение непрерывно зависит от входных данных  h и  h , и эта зависимость равномерна относительно шага сетки, т.е. для любого   0 существует такое  ( ) , не зависящее от шага h (по крайней мере, для достаточно малых h), что y (1)  y ( 2 )   , если  h(1)   h( 2)   ,  h(1)   h( 2)   . В случае нескольких независимых переменных рассматривают условную и безусловную устойчивость. Если сформулированные выше условия выполняются при определенном соотношении шагов по различным переменным, то устойчивость называется условной, в противном случае, когда соотношение между шагами может быть произвольным, устойчивость называется безусловной. При этом непрерывную зависимость решения от  h называют устойчивостью по правой части, а непрерывную зависимость от  h - устойчивостью по дополнительным условиям (начальным и граничным). Устойчивость разностных схем может быть исследована несколькими методами: разделения переменных, энергетических неравенств, операторных неравенств, на основе принципа максимума и др. Рассмотрим вначале принцип максимума. Перепишем разностные схемы (14) в виде a k k y n k   b p y n p   n , p (23) где суммирование выполняется по узлам шаблона около n-го узла. Пронумеруем узлы так, чтобы a0  max a k . k Принцип максимума, дающий достаточное условие устойчивости явных и неявных двухслойных разностных схем, формулируется следующим образом: 1) схема равномерно устойчива по начальным данным, если (1  C ) a0   a k   b p , C  const  0 . k 0 (24) p 2) схема устойчива по правой части, если справедливо соотношение (24) и имеет место неравенство a0   ak  k 0 D  , D  const  0 (25) Сформулированные условия не являются необходимыми для устойчивости схем, т.е. несоблюдение (24), (25) не обязательно ведет к их неустойчивости. Данным методом можно доказать устойчивость схем точности O( ) , в других случаях используют иные методы. В качестве примера рассмотрим неявную схему (7). Представив ее в виде (23), получим выражения для коэффициентов a0  2a 1 a 1  , a 1  a1  2 , b0  2   h h для 1  n  N-1. В случае граничных условий первого рода (10) a 0  1 ,  0  0 для n=0 и n=N. Все остальные коэффициенты равны нулю. Видно, что условие (25) выполнено во внутренних узлах при произвольных соотношениях шагов по переменным x и t, а условие (24) справедливо во всех узлах сетки. Таким образом, рассматриваемая неявная разностная схема безусловно устойчива по начальным данным, правой части и краевым условиям. Принцип максимума позволяет доказывать устойчивость в чебышевской (локальной) норме. Метод разделения переменных применяют для исследования устойчивости схем в гильбертовой (среднеквадратичной) норме. Запишем разностную схему в канонической форме B y  y где  A  Ay   , и B- (26) разностные операторы, действующие на функцию по пространственным переменным. Например, для явной схемы (8) очевидно, что B  E , Ay  a y n1  2 y n  y n1 . h2 (27) При точной правой части погрешность решения удовлетворяет уравнению Bz  ( A  B) z  0 . (28) Частное решение (28) будем искать методом разделения переменных: z ( x n , t m )   qm e iqxn / l , q=0,  1,  2, … (29) Здесь  q - множитель роста q-ой гармоники при переходе с одного слоя на следующий слой, так что z   q z . Подставляя (29) в (28), придем к уравнению  q Be iqx / l  ( A  B)e iqx / l  0 . (30) Если в схеме (26) коэффициенты постоянны, а сетка равномерна, то уравнение (30) после сокращения множителя e iqx / l не будет зависеть от индекса n, т.е. от координаты x, соответственно  q не будет зависеть от x. Признак устойчивости формулируется следующим образом. Схема (26) с постоянными коэффициентами устойчива по начальным данным, если для всех гармоник с индексом q выполняется неравенство  й  1  С , C  const Константа C не должна быть большой, поэтому обычно принимают C=0. (31) Признак неустойчивости заключается в следующем: если хотя бы для одной гармоники q величину  q нельзя мажорировать величиной 1  C , то схема (26) неустойчива. В качестве примера проверим устойчивость явной схемы (8). Зафиксируем правую часть в разностной схеме (26), представим погрешность решения в виде (29) с учетом того, что z   q z , и, подставив ее в (8), получим  (  q  1)e iqxn / l      m q iq ( xn  h ) / l   2e iqxn / l  e iq ( xn  h ) / l   a qm  e   h2      0 ,  или (  q  1)  a iqh / l (e  2  e iqh / l )  0 . 2 h Заменяя по формуле Эйлера e iqh / l  eiqh / l  2 cosqh / l , получаем окончательно (  q  1)  4 a qh sin 2  0. 2 2l h Отсюда q  1  4 a qh sin 2 2 2l h (32) Теперь согласно условию (31) (  q  1 ) получим из (32) критерий устойчивости рассматриваемой разностной схемы, учитывая, что sin 2 qh 2l 0 4a 2. h2 Таким образом, явная схема устойчива при определенном соотношении между шагами по координатам x и t:  h2 , 2a (33) т.е. явная схема условно устойчива. Отметим, что из признака (31) и дополнительного условия (25) следует устойчивость разностной схемы по правой части в гильбертовой норме. Если дифференциальное уравнение имеет переменные коэффициенты или используется неравномерная сетка, то применяют прием «замораживания коэффициентов». В этом случае коэффициенты считают постоянными и равными их значениям в некотором фиксированном узле n. Схема считается устойчивой, если при любых n и q оказывается справедливым неравенство (31). Применим метод разделения переменных для исследования устойчивости разностной схемы (9). Вначале проверим устойчивость по начальным данным. Положим   0 и представим погрешность решения в виде (29), при этом z   z . Сделав данную n q подстановку в (9), получим  (  q  1)e iqxn / l     e  a(1   )  iq ( xn  h ) / l   2e iqxn / l  e iq ( xn  h ) / l   a q  e   h2   iq ( xn  h ) / l  2e iqxn / l h e iq ( xn  h ) / l 2    0      , или, проведя сокращение на e iqxn / l , перейдем к уравнению (  q  1)   e  iqh / l  2  e iqh / l  a q  h2   e iqh / l  2  e iqh / l  a(1   ) h2     0      . Используя формулу Эйлера, получаем окончательно 4a qh 4a (1   ) 2 qh  q sin 2  sin  0. 2 2l 2l h h2 q 1  Отсюда q  4a (1   ) 2 qh sin 2l . h2 4a  qh 1  2 sin 2 2l h 1 (34) Из этого выражения следует, что при любом   0 множитель роста гармоники  q  1. Осталось выяснить условие выполнения соотношения  q  1. Элементарными преобразованиями легко показать, что данное обстоятельство реализуется в случае, если 4 a (1  2 )  2 h2 или  1 h2  2 4a (35) Дополнительное условие устойчивости по правой части (25), как это следует из (16), выполняется для всех  и h. Таким образом, разностная схема (9) устойчива по правой части, если выполнено условие (35) устойчивости по начальным данным. Для чисто неявной разностной схемы   1и из (35) следует, что данные схемы устойчивы при любом соотношении шагов по независимым переменным, т. е. они безусловно устойчивы. Этот результат был получен выше на основе принципа максимума. Для явной разностной схемы   0 , и устойчивость согласно (35) обеспечивается при условии   h2 , о чем было сказано ранее. 2a 17.2. Сходимость Дадим определение сходимости. Разностное решение y(x) сходится к решению задачи (17), (18), если y( x )  u( x )  0 при h  0 ; разностное решение имеет порядок точности p, если y( x )  u( x)  O( h p ) при h  0 . В теории разностных схем большое значение имеет теорема, которую часто в кратком виде формулируют следующим образом: «Из аппроксимации и устойчивости следует сходимость». Сформулируем и докажем данную теорему. Теорема. Если решение задачи (17), (18) существует, разностная схема (19), (20) корректна (т. е. ее решение существует, единственно и схема устойчива) и аппроксимирует задачу на данном решении, то разностное решение сходится к точному. Действительно, из определения невязки (21), (22) следует Ah u   h   , (36) Bh u   h   . (37) Перепишем разностную схему (19), (20) Ah y  h , Bh y   h . Сравнивая (36), (37) с выписанной разностной схемой, видим, что система (36), (37) представляет собой не что иное, как разностную схему, правые части которой изменены на величину невязки. Из устойчивости разностной схемы следует, что   0  ( ) , такое, что y  u   , если    ( ),    ( ) . Но разностная схема аппроксимирует задачу, значит   0 h0 ( ) , такой, что невязки    ,    при h  h0 ( ) . Таким образом,   0 h0 ( ( )) , такой что y  u   при h  h0 ( ) , что и доказывает сходимость. Сделаем замечание относительно точности разностного решения. Краткая формулировка соответствующей теоремы имеет следующий вид: «Для линейных разностных схем порядок точности не ниже порядка аппроксимации». В развернутом варианте теорема формулируется так: Если условия вышеприведенной теоремы о сходимости выполнены, операторы Ah , Bh линейные, а порядок аппроксимации равен p , то сходимость имеет порядок не ниже p . ВОПРОСЫ ДЛЯ САМОКОНТРОЛЯ 1. Какие дифференциальные уравнения называются уравнениями в частных производных? 2. Приведите классификацию уравнений второго порядка, линейных относительно старших производных. 3. Какие существуют постановки задач для уравнений указанного типа? 4. Какие типы граничных условий можно сформулировать? 5.Дайте определения многомерной разностной сетки, узлов, слоев, направлений, шаблона, разностной схемы. 6. Постройте явную и неявную разностную схему одномерного уравнения параболического типа на трехточечном и шеститочечном шаблонах. 7. Дайте определение аппроксимации, устойчивости и сходимости разностных схем. 8. Оцените порядок аппроксимации разностной схемы, получив выражение для невязки в случае шеститочечной разностной схемы для уравнения параболического типа. 9. Проверьте аппроксимацию граничных условий третьего рода в простейшем варианте построения разностного аналога производной на основе односторонней разности. 10. Оцените устойчивость явных и неявных разностных схем для параболического уравнения, используя принцип максимума и метод разделения переменных. 11. Докажите теорему о сходимости приближенного разностного решения к точному решению исходной дифференциальной задачи. 12. Как соотносятся порядок аппроксимации и порядок точности разностных схем.