Вычислительная математика лекции. Конспект лекций. Санкт петербург 2011 2 Оглавление
 Скачать 3.86 Mb.
|
|
15.2.8.1. Простейший вариант метода Ньютона. Пусть в заданной области ||x-x*||<ϭ выполнены необходимое и достаточное условия существования минимума функции f(x). Тогда задача вычисления x*=arg min x f(x) сводится к нахождении корней системы нелинейных уравнений g(x)=0, где g(x)= grad(f(x)). Для решения задачи можно использовать метод Ньютона g(x n )+G(x n )(x n+1 -x n )=0. Здесь x n n – ое приближение, G(x n ) матрица Гессе. Таким образом, чтобы попасть из точки x n в точку x n+1 , нужно переместится вдоль вектора p n =x n+1 -x n , который определяется из системы линейных алгебраических уравнений G(x n )p n =-g(x n ). Вектор p n называют ньютоновским направлением, а метод спуска x n+1 =x n +α n p n – методом Ньютона. Ньютоновское направление действительно является направлением спуска. В самом деле, p n =-G -1 (x n )g(x n ). Матрица G -1 (x n ) положительно определена (это следует из положительной определенности матрицы G(x n )). Поэтому (g(x n ),p n )=-(G -1 (x n )g(x n ),g(x n )<0. При выборе α n =1, метод в точности совпадает с методом Ньютона решения системы нелинейных уравнений g(x)=0. Как было выяснено в соответствующем разделе, методу в таком варианте присущ существенный недостаток, связанный с необходимостью тщательного подбора начального приближения. Его достоинство – квадратичная скорость сходимости. 193 Метод Ньютона может быть интерпретирован как последовательная квадратичная аппроксимация. Действительно, при выполнении условий непрерывности любую функцию f(x) можно аппроксимировать в окрестности точки x 0 функцией 0 0 0 0 0 0 1 ( ) ( ) ( ) ( ) ( ) ( )( ) 2 T x f x x x f x x x G x x x , где G(x 0 ) – матрица Гессе, вычисленная в точке x 0 Разумной аппроксимацией минимума функции f(x) может быть минимум функции φ(x). Если последний находится в точке x m , то 0 0 0 '( ) 0 или ( ) ( )( ) 0, m m x f x G x x x откуда 1 0 0 0 1 0 0 0 - ( ) ( ), - ( ) ( ). m m x x G x f x или x x G x g x Полагая x 0 =x n , x m =x n+1 , получаем формулу метода Ньютона. 1 1 - ( ) ( ). n n n n x x G x g x 15.2.8.2. Метод Ньютона с дроблением шага. Позволяет облегчить проблему выбора начального приближения, однако не обеспечивает достаточную скорость вычисления. Аналогично методу наискорейшего спуска очередной шаг α n выбирается из условия 1 ( ) ( ) ( , ), 0< <1. n n n n f x f x g p Теорема о сходимости. Пусть выполнены следующие условия: 1.f(x) ограничена снизу. 2 Матрица Гессе удовлетворяет условиям 2 2 || || ( ( ) , ) || || ; , 0 y G x y y R y R 3. 1 1 ( ) ( ) ( , ), 0< <1, где ; n n n n n n n f x f x g p p G g x n+1 = x n +α n p n . 194 Тогда ||g n || → 0 при n →∞ . Целевая функция f(x) монотонно уменьшается с ростом n. В процессе доказательства теоремы формулируются требования к выбору коэффициента 2 2 (1 ). Алгоритм метода, учитывающий эти условия, выглядит следующим образом. 1. Выбираем произвольную точку x 0 – начальное приближение. Полагаем x n =x 0 2. Выбираем произвольное значение α. 3. Определяем новую точку x n+1 =x n +αp n =x n -αG -1 (x n )g(x n ). . 4. Вычисляем f(x n+1 )=f(x n +αp n ) и f(x n ). 5. Если f(x n+1 )-f(x n ) γα(g(x n ), p n ), 0<γ<1, то значение α берем в качестве искомого α n и переходим к шагу 6, иначе уменьшаем α ( путем умножения на множитель γ, 0<γ<1) и переходим к шагу 4. 6. Вычисляем новое приближение x n+1 =x n -αG -1 (x n )g(x n ). . 7. Проверяем выполнение условия стационарности ||g n+1 ||= (g n+1 ,g n+1 ) ε – заданная погрешность. Если условие выполнено переходим к шагу 8, иначе к шагу 2, положив x n =x n+1 8. Останов. При использовании метода Ньютона проблемы возникают в случае плохой обусловленности матрицы Гессе. Эта ситуация была названа ранее проблемой овражности. 15.2.8.3. Понятие о квазиньютоновских методах. В методах, которые называются квазиньютоновскими, также как и в методах Ньютона используется квадратичная аппроксимация целевой функции. Однако, здесь не производится вычисление матрицы Гессе. Информация о кривизнах функции здесь накапливается постепенно в результате наблюдений за изменением градиента. 195 Направление спуска в квазиньютоновских методах определяется как решение системы уравнений ( ) ( ) ( ) , n n n B p g в которой ( ) n B - текущее приближение к матрице Гессе. Начальное приближение обычно полагают равным единичной матрице. В таком случае p (0) = -g (0) и первая итерация совпадает с одним шагом метода наискорейшего спуска. После того, как найдено приближение x (n+1) вычисляют матрицу B (n+1) = B (n) + ΔB (n) , где ΔB (n) некоторая поправочная матрица. В пределе матрица ( ) n B становится матрицей Гессе. Среди наиболее популярных квазиньютоновских алгоритмов известен BFGS – алгоритм (Broyden, Fletcher, Goldfarb, Shanno), в котором аппроксимация ( 1) n B производится итерационно по формуле ( ) ( ) ( ) ( ) ( ) ( ) ( 1) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) , ( ) ( ) n n T n n n T n n n n T n n T n n q q B s s B B B q s s B s где s (n) = x (n+1) – x (n) , q (n) = g (n+1) - g (n) В некоторых вариантах квазиньютоновских методов аппроксимируется матрица ( ) n H , обратная ( ) n B (метод Давидона - Флетчера – Пауэла). Рассматриваемые методы при определенных условиях сходятся сверхлинейно, однако весьма чувствительны к ошибкам расчета градиента. Квазиньютоновские методы непрерывно совершенствуются и к настоящему времени стали одними из наиболее широко применяемых на практике. В этих методах сочетаются как идеи метода Ньютона – Рафсона, так и свойство сопряженных направлений. Подробнее с материалом можно ознакомиться в следующих источниках: 1. Гилл Ф., Мюррей У., Райт М. Практическая оптимизация. – М.; Мир, 1985. 196 2. Банди Б. Методы оптимизации. Вводный курс. – М. : Радио и связь. 198 15.2.9. Метод Левенберга – Марквардта. Рассматриваемый метод является модификацией метода Ньютона, направленной на устранение недостатка, связанного с необходимостью тщательного подбора начального приближения, и частично проблемы овражности Направление поиска Левенберга — Марквардта определяется из системы: (G(x n )+λ n E) p n =-g(x n ). где λ n — некоторая неотрицательная константа, своя для каждого шага, E — единичная матрица, x n+1 =x n + p n Нетрудно заметить, что при λ n = 0 алгоритм вырождается в метод Ньютона. При достаточно большом λ n метод незначительно отличается от метода наискорейшего спуска, сходимость которого гарантирована при произвольном выборе начального приближения. Таким образом, при правильном подборе параметра λ n можно добиться сочетания достоинств методов наискорейшего спуска и Ньютона . Величину λ n можно увеличивать, пока выполняется условие f(x n+1 ) < f(x n ) . Однако при чрезмерно большом λ n проявляются все недостатки метода градиентного спуска, связанные с проблемой овражности. Дополнительное решение этой проблемы предложил Марквардт. Он заметил, что если заменить единичную матрицу на диагональ матрицы Гессе, то можно достичь увеличения шага вдоль пологих участков и уменьшения вдоль крутых спусков: (G(x n )+λ n diag(G(x n )) p n =-g(x n ). 197 16.Методы решения краевой задачи для линейных обыкновенных дифференциальных уравнений. 16. 1. Постановка задачи. Ранее уже рассматривались методы решения обыкновенных дифференциальных уравнений – методы решения задачи Коши: (t) = f(x,t), x(t 0 ) = x 0 , t [t 0 ,t k ] , x(t) = ? Краевая задача отличается от данной заданием краевых вместо начальных условий, то есть условий в более чем одной точке внутри интервала интегрирования: (t) = f(x,t), t [a,b] , α T x(a) = γ a , β T x(b) = γ b , где α, β – заданные числовые векторы одной размерности с x, γ a , γ b – заданные числа. Граничные условия могут быть заданы и в более сложной форме. В отличие от имеющей единственное решение задачи Коши краевая задача может иметь или одно решение, или бесконечное множество решений, или, наконец, может совсем не иметь решений. Рассмотрим наиболее часто встречающуюся на практике краевую задачу второго порядка: Lx ≡ (t) + p(t) (t) + q(t) x(t) = f(t), (1) и краевые условия αx ≡ α 0 x(a) + α 1 (a) = γ 0 , (2 а ) βx ≡ β x(b) + β 1 (b) = γ 1 , (2 в ) где p(t), q(t), f(t) C [a,b] – заданные функции, α i , β i , γ i – заданные числа. Теорема. Для того, чтобы существовало единственное решение неоднородной краевой задачи (1), (2) необходимо и достаточно, чтобы 198 однородная краевая задача Lx ≡ 0; αx ≡ , βx ≡ , t [a,b] имела только тривиальное решение x(t) ≡ Однако и это условие не всегда удается проверить. На практике существует большое число методов численного решения краевой задачи. Два из них рассмотрены ниже. 16. 2. Метод стрельбы (баллистический метод). Рассмотрим метод на примере следующей задачи: (t) = f(x,t), t [a,b] , x R 2 , или 1 (t) = f 1 (t,x 1 , x 2 ) (3) 2 (t) = f 2 (t,x 1 , x 2 ) с граничными условиями φ(x 1 (a), x 2 (a)) = 0, (4 a ) ψ(x 1 (b), x 2 (b)) = 0. (4 b ) Выберем произвольно x 1 (a) = η и найдем x 2 (a), решив уравнение (4 a ): x 2 (a) = ζ(η). После этого решим задачу Коши для системы (3), используя начальные условия x 1 (a) = η и x 2 (a) = ζ(η), t [a,b]. На правой границе найдем x 1 (b,η) и x 2 (b,η). При этом почти наверняка условие (4 b ) не будет выполнено ψ(x 1 (b,η), x 2 (b,η)) = ψ (η) ≠ 0. Очевидно, нужно изменить левое краевое условие так, чтобы ψ (η) = 0. Тем самым решение краевой задачи (3), (4) свелось к решению одного алгебраического уравнения ψ (η) = 0. (5) Сложность в том, что это уравнение задано не аналитически, а алгоритмически. Простейший метод решения уравнения (5) состоит в следующем. Делают пробные "выстрелы" – расчеты с наудачу выбранными значениями η i до тех пор, пока ψ (η i ) и ψ (η i+1 ) окажутся разных знаков. Далее, используя метод половинного деления на отрезке [η i , η i+1 ], продолжаем процесс решения (5) до получения требуемй точности |ψ (η n ) | < ε. 199 Рассматриваемый метод весьма трудоемок, так как требует многократного интегрирования системы (3). 16.3. Разностный метод (метод конечных разностей). Рассмотрим краевую задачу: Lx ≡ (t) + p(t) (t) - q(t) x(t) = f(t), (1) t [0,1] , x(0) = x 0 , x(1) = x 1 , (2) где p(t), q(t), f(t) C 2 [0,1] – заданные функции, q(t) ≥0 , x 0 , x 1 – заданные числа. Теорема. При выполнении указанных условий задача (1), (2) имеет единственное решение x(t) C 4 [0,1] Решение задачи предполагает использование сетки сеточных функций. На отрезке [0,1] задается конечное множество точек (узлов) с постоянным шагом {t k } N k=0 , где t k =kh, h=1/N, N ≥2. Значения сеточной функции в узлах будем обозначать x k Аппроксимируем значения производных соответствующими разностными операторами 1 1 1 1 2 '( ) , 2 2 ''( ) k k k k k k k x x x t h x x x x t h Тогда вместо дифференциальной задачи (1), (2) получаем разностную 0 1 1 1 1 1 0 1 2 2 , k=1, 1; ; x 2 k k k k k k k k k x x x x x p q x f N x x x h h Или b k x k-1 + c k x k + d k x k+1 , (2) где 2 2 2 0 1 0 1 1 2 1 ; c ; d ; 2 2 при условии k=1, 1; ; x k k k k k k p p b q h h h h h N x x x Разностное уравнение (2) при всех значениях k эквивалентно системе линейных алгебраических уравнений 200 b 1 x 0 + c 1 x 1 +d 1 x 2 = f 1 , b 2 x 1 + c 2 x 2 +d 2 x 3 = f 2 , ……………… b N-1 x N-2 + c N-1 x N-1 +d N-1 x N = f N-1 Если дополнить эти уравнения сверху соотношением x 0 = x 0 и снизу x N = x 1 , то получим линейную алгебраическую систему Ax = r , где x = (x 0 ,x 1 ,…,x N ) T , r = (x 0 ,f 1 ,f 2 ,…,f N-1 , x 1 ) T и А – трёхдиагональная матрица 1 1 1 2 2 2 1 1 1 1 0 0 1 N N N b c d b c d A b c d Если выполняется условие |c k | ≥|b k | + |d k |, то говорят, что матрица А является матрицей с доминирующей главной диагональю. Указанное условие гарантирет существование и единственность решения. Оно всегда выполняется, если |p k |h 2. Наиболее экономично решать задачу с трехдиагональной матрицей методом прогонки. 17. Сингулярное разложение и метод наименьших квадратов. 17. 1. Решение переопределенной системы уравнений. Вернемся к задаче среднеквадратичного приближения сеточной функции. Заданную табличную функцию {f j , x j } размера m (j=0,1,…,m-1) требуется заменить обобщенным полиномом степени (n-1), имеющим набор базисных функций φ i (x), i=0,1,…,n-1. 201 1 1 0 ( ) ( ) n n i i i P x x Обычно m>n. В такой постановке задача сводится к попытке решения системы линейных алгебраических уравнений Aα=f, A R m×n , α R n , f R m 0 0 1 0 0 0 0 1 1 1 1 1 1 0 ( ) ( ) , = , f= ( ) ( ) ( ) , j=0,1,...,m-1 (1) n m n m n m n i j i j i x x f A x x f x f Указанная система может иметь решение и при том решение единственное, если при m>n вектор f лежит в пространстве столбцов матрицы А и матрица является матрицей полного ранга, т. е. r=min(m,n)=n. Для квадратной матрицы выполнение последнего условия означает её невырожденность. Когда r 2 2 0 ( ) ( ) || || m j i i j j F x f x A f Таким образом, следует минимизировать норму невязки ||Aα-f|| 2 или длину вектора Aα-f. Длина этого вектора есть расстояние от точки f до пространства столбцов. Минимальной длиной обладает перпендикуляр, опущенный из точки f на пространство столбцов. Таким образом, вектор невязки, удовлетворяющий условиям среднеквадратичного приближения, должен быть ортогонален к этому пространству. 202 Столбец 1 A α Столбец 2 f f-A α На рисунке изображен случай n=2, m=3. Предполагается, что столбцы матрицы А линейно независимы, т. е. А матрица полного ранга. 0 0 1 0 0 1 1 1 0 2 1 2 ( ) ( ) 1 ( ) , столбец2= ( ) ( ) ( ) x x столбец x x x x Проекция p вектора f на пространство столбцов это вектор, лежащий в пространстве столбцов. Решение системы Aα=p существует, единственно и дает искомый вектор α, удовлетворяющий критерию наименьших квадратов. Ортогональность вектора невязки к пространству столбцов состоит в следующем. Каждый вектор в пространстве столбцов матрицы А является линейной комбинацией столбцов с некоторыми коэффициентами y 0 ,…y n-1 , т. е. это вектор вида Ay=A (0) y 0 +…+A (i) y i +…+A (n-1) y n-1 , где A (i) - i – ый столбец матрицы А, y i - i – ый компонент вектора y. Вектор Ay и вектор невязки должны быть взаимно перпендикулярны при любом y: (Ay,Aα-f)=(Ay) T (Aα-f)=0, или y T (A T Aα-A T f)=0. Так как y произволен, то вектор в скобках должен быть нулевым. Решение по методу наименьших квадратов для несовместной системы Aα=f удовлетворяет соотношению A T Aα = A T f. Оно известно под названием нормальных уравнений. Проекция точки f на пространство столбцов имеет вид 203 p= Aα=A(A T A) -1 A T f. Матрица (A T A) -1 существует тогда и только тогда, когда матрица А является матрицей полного ранга. Тот же самый результат можно получить аналитически следующим образом. Нужно минимизировать квадрат длины вектора невязки (r,r)=(Aα-f,Aα-f) →min. (r,r)=(Aα-f,Aα-f)=(Aα,Aα)-2(Aα,f)+(f,f)=(A T Aα,α)-2(α,A T f)+(f,f). Матрица A T A симметрическая размера n. Из условия минимума: ( , ) 2 2 0 T T r r A A A f следует соотношение A T Aα = A T f, называемое системой нормальных уравнений. Если столбцы матрицы А линейно независимы, то матрица A T A R n×n обратима и единственное решение определяется формулой α =(A T A) -1 A T f. Линейная зависимость столбцов матрицы А возникает, если базисные функции не обладают свойством независимости. В этом случае критерий наименьших квадратов необязательно определяет единственный набор коэффициентов. Действительно, линейная зависимость базисных функций в заданных точках означает существование коэффициентов γ i , среди которых есть ненулевые, таких, что 1 0 ( ) 0, j=0,1,...,m 1, или, что тоже самое, A =0. n i i j i x Тогда, как следует из (1), коэффициентам α i можно прибавить любое кратное коэффициентов γ i , не изменяя длину вектора невязки. Важной и часто забываемой задачей при выравнивании данных методом наименьших квадратов с произвольными базисными функциями является обнаружение такой зависимости и принятие надлежащих мер. Таким образом, задача наименьших квадратов все еще остается не до конца решенной. Как быть, если матрица A T A окажется необратимой? 204 Дополнительно нужно учитывать, что при n>5 матрица A T A становится настолько плохо обусловленной, что её обращение становится практически невозможным, даже если условия независимости базисных функций соблюдаются. Универсальным способом решения задачи наименьших квадратов является использование сингулярного разложения. |