Вычислительная математика лекции. Конспект лекций. Санкт петербург 2011 2 Оглавление
 Скачать 3.86 Mb.
|
|
Р е к у р р е н т н а я а п п р о к с и м а ц и я м а т р и ч н о й э к с п о н е н т ы. Решение системы разностных уравнений требует вычисления матриц S = e AH и G= 0 H A e Bd , которые не зависят от величины t k Их вычисление может быть проведено с использованием рекуррентных формул. Пусть h=H/2 N , где N целое число. Вычислив exp(Ah), можно найти exp(AH), применяя последовательно формулы exp(2Ah)= exp(Ah)exp(Ah) ………………………………. exp(AH)= exp(2 N Ah)=exp(2 N-1 Ah)exp(2 N-1 Ah) Таким образом, получаем рекуррентное соотношение φ k+1 =φ k 2 , где φ k =exp(2 к Ah), k=0,1,…,N-1, причем φ N =exp(AH)=S. Построим теперь рекуррентное соотношение для вычисления матрицы G= 0 H A e Bd Докажем следующее соотношение 2 0 0 exp( ) exp( ) ( exp( )) h h A d A d E Ah . С этой целью преобразуем правую часть 0 0 0 2 2 0 0 2 0 exp( ) ( exp( )) exp( ) exp( ) exp( ) exp( ) exp( ) ; exp( ) exp( ) h h h h h h h h h h h A d E Ah A d A Ah d A d A d A d A Ah d A d Таким образом, g 1 =(E+φ 0 )g 0 , где 98 1 2 0 1 0 0 1 2 2 1 0 0 exp( ) ; exp( ) ; ( ); k=0,1,...,N-1; exp( ) ; exp( ) ; k k h h k k k h h k k N A d B g A d B g g g E A d B g A d B g g G О в ы б о р е ш а г а h. Шаг h должен выбираться достаточно малым для достижения приемлемой точности аппроксимации при разумном числе l учитываемых слагаемых в разложении экспоненты. Обозначим φ 0 сумму l учитываемых слагаемых 1 0 1 2 2 1 1 1 || ||/ ( ) ( ) exp( ) ! ( 1)! (|| ||) (|| || ) || || ! ( 1)! (|| ||) || || (|| || ) ( || ||) || || (1 ...) ,( 1). ! ! l l R l l l l w h A l Ah Ah Ah l l Ah A h R l l Ah A h A h h A h A l l l l l Таким образом, погрешность аппроксимации ( || ||) ! l h A l . При заданной величине δ должно быть h< ! || || l l A Шаг дискретности H=2 N h выбирается из условия наблюдения процесса, описываемого дифференциальным уравнением. При этом необходимая точность обеспечивается малым значением начального шага h. Использование рекуррентных соотношений дает выигрыш по количеству операций в ≈2 N раза в сравнении с применением малого шага h на всем диапазоне решения. 99 9.12. Матричные нормы. Определение матричной нормы эквивалентно определению векторной нормы. Именно, функция f, отображающая элементы множества R m×n в элементы множества вещественных чисел R, является матричной нормой, если f(A)=||A|| ≥0 A R m×n (||A||=0 ⟺A=0), ||A+B|| ||A||+||B|| A, B R m×n , ||αA||=|α| ||A|| α R, A R m×n Чаще всего используются нормы Фробениуса (евклидова норма) 2 1 1 || || m n F ij i j A a и p – нормы ||A|| p = 0 || || sup || || p x p Ax x О p – нормах говорят, что они подчинены соответствующим векторным нормам. Для таких норм выполняется условие ||Ax|| p ||A|| p ||x|| p Введенные матричные нормы удовлетворяют мультипликативному условию ||AB|| ||A|| ||B||. При выполнении этого условия говорят, что нормы матриц A и B взаимно согласованы. Некоторые свойства матричных норм. Для A R m×n имеем: ||A|| 2 ||A|| F 2 || || , n A 2 , , | | || || | |, max max ij ij i j i j a A mn a 100 ||A|| 1 = 1 1 | | max m ij j n i a , ||A|| ∞ = 1 1 | |, max n ij i m i a 2 1 || || || || || || , A A m A n 1 2 1 1 || || || || || || , A A n A m 2 1 || || ( ). max T j j n A A A 9.13. Устойчивость решений обыкновенных дифференциальных и разностных уравнений. 9.13.1. Основные понятия теории устойчивости. Факт устойчивости устанавливается по результатам анализа решений исходной и возмущенной систем уравнений. Пусть исходная и возмущенная системы имеют вид: 0 0 0 0 . система ( ( ), ( )), нач. усл. x(t ) * . система ( * ( ), ( )), нач. усл. x*(t ) * dx исх f x t u t x dt dx возм f x t u t x dt Обозначим x(t)-x*(t)=v(t). О п р е д е л е н и е. Решение x(t) исходной системы устойчиво по Ляпунову, если для любого ε > 0 существует δ > 0, такое что || v(t) || ε для всех t > t*, если || x 0 – x 0 * || δ. 101 t* -ε ε δ ν (t) t Если решение x(t) - устойчиво по Ляпунову - норма || v(t) || → 0 при t→ ∞, то это решение асимптотически устойчиво. Аналогичная ситуация имеет место при анализе устойчивости решения систем разностных уравнений Д о п о л н и т е л ь н ы е з а м е ч а н и я. 1.Устойчивость по Ляпунову эквивалентна выполнению условия 0 0 * || ( ) *( ) || || * || max t t x t x t K x x Скаляр K играет роль числа обусловленности. Это условие получается, если положить ε=Кδ. 2.Возмущенная система может возникать при внесении погрешностей как в начальные условия, так и в правую часть. Можно показать, что в этом случае для устойчивой по Ляпунову системы выполняется условие 102 0 0 0 * || ( ) *( ) || [|| * || || ( ) ||] max max t t t x t x t K x x . Здесь ψ(t) малая вектор-функция для дифференциальной системы или малая сеточная вектор-функция для разностной системы, добавленная в виде слагаемого к правой части. 9.13.2..Анализ устойчивости решения систем линейных обыкновенных дифференциальных уравнений с постоянными коэффициентами. Исходное уравнение имеет вид ( ); A , B ; x , f n n n m n m dx Ax Bf t dt x(t 0 ) = x 0 , t [t 0 , T]. Полагая как и прежде x(t)-x*(t)=v(t), x(t 0 )- x*(t 0 )=v(t 0 ) = v 0 ,получаем 0 0 , v(t ) dv Av v dt Здесь x*(t) решение возмущенной системы, полученное из исходной при новых начальных условиях. Таким образом, для установления факта устойчивости нужно проанализировать решение однородной системы дифференциальных уравнений. Это решение вычислено ранее 0 ( ) 0 v(t)=e v(t ). A t t Для вычисления матричной экпоненты воспользуемся полиномом Лагранжа – Сильвестра. Пусть минимальный аннулирующий многочлен матрицы А имеет вид ψ(λ) = (λ-λ 1 ) m 1 (λ-λ 2 ) m 2 …(λ-λ s ) m s , m 1 + m 2 +…+ m s =m.В данном случае f(λ) = 0 ( ) t t e Тогда 103 0 1 0 0 0 1 0 0 ( ) ( 1) 1 2 1 ( ) ( ) ( ) ( 1) 1 2 1 ( ) 1 ( ) 1 2 0 0 1 [ ( ) '( ) ... ( )] [ ( ) ' ... ( ) ] [ ( ) ... ( ) ] k k k k k k k k k k s A t t m k k k k km k k s t t t t t t m k k km k s A t t m t t k k km k e Z f Z f Z f Z e Z e Z e e Z Z t t Z t t e Здесь Z ki матрицы с постоянными коэффициентами. Рассмотрим три случая: 1.Reλ k 0 (k = 1,2,…,s), причем для тех λ k , для которых Reλ k =0, m k = 1 (т. е. чисто мнимые характеристические числа являются простыми корнями минимального многочлена). 2. Reλ k <0 (k = 1,2,…,s). 3.Хотя бы для одного k имеем Reλ k >0 либо Reλ k =0, но m k > 1.Тогда из формулы 0 0 ( ) 1 ( ) 1 2 0 0 1 [ ( ) ... ( ) ] k k k s A t t m t t k k km k e Z Z t t Z t t e следует, что в случае (1) элементы матрицы Q(t)= 0 ( ) A t t e ограничены на интервале [t 0 , ∞), в случае (2) Q(t)→0 при t→∞, а в случае (3) элементы матрицы Q(t) не ограничены на интервале (t 0 , ∞). Таким образом, в случаях (1) и (2) решение устойчиво по Ляпунову причем в случае (2) устойчиво асимптотически. В третьем случае неустойчиво. 9.13.3.Анализ устойчивости решения систем линейных разностных уравнений с постоянными коэффициентами. Исходное уравнение имеет вид 1 max ; S , G ; x(t )=x , f(t )=f [0, ]. n n n m n m k k k k k k k x Sx Gf k k Начальные условия: x(t 0 ) = x 0 . Если возмущать начальные условия, обозначив, как и прежде, x k -x* k =v k , x 0 -x* 0 = v 0 ,получаем 1 k k v Sv 104 Таким образом, для установления факта устойчивости нужно проанализировать решение однородной системы разностных уравнений. Это решение получено ранее v k = S k v 0 Пусть минимальный аннулирующий многочлен матрицы S имеет вид ψ(λ) = (λ-λ 1 ) m 1 (λ-λ 2 ) m 2 …(λ-λ s ) m s , m 1 + m 2 +…+ m s =m. В данном случае f(λ) = λ k Тогда 1 1 2 1 [ ( ( 1)...( 1)) ]. i i s k m k k k i i i i im i i i S Z Z k Z k k k m Рассмотрим три случая: 1.| λ i | 1 (i = 1,2,…,s), причем для тех λ i , для которых | λ i | =1, m i = 1 (т. е. соответствующие характеристические числа являются простыми корнями минимального многочлена). 2. | λ i | <1 (i = 1,2,…,s). 3.Хотя бы для одного i имеем | λ i | >1 либо | λ i | =1, но m i > 1.Тогда из формулы 1 1 2 1 [ ( ( 1)...( 1)) ] i i s k m k k k i i i i im i i i S Z Z k Z k k k m следует, что в случае (1) матрица S k ограничена на интервале k [0,∞), в случае (2) S k →0 при k→∞, а в случае (3) матрица S k не ограничена на интервале k [0,∞). Таким образом, в случаях (1) и (2) решение устойчиво по Ляпунову, причем в случае (2) устойчиво асимптотически. В третьем случае неустойчиво. 105 10. Численные методы решения обыкновенных дифференциальных уравнений 10. 1. Постановка задачи. Рассмотрим численные методы решения известной задачи Коши: 0 0 0 ( , ), t [t , ], x(t ) , x n k dx f t x t x dt Прежде чем начинать вычисления, необходимо установить их принципиальную возможность, т. е. факт корректности задачи. Иными словами нужно установить 1) существование решения, 2) единственность решения, 3)непрерывную зависимость погрешности решения от погрешности исходных данных (в данном случае начальные условия и правая часть). Ответ на эти вопросы дает теорема Коши существования и единственности. Если f(t,x) непрерывна в некоторой области вокруг точки (t 0 ,x 0 ), тогда существует по крайней мере одно решение уравнения x'=f(t,x), принимающее при t=t 0 значение x 0 , определенное и непрерывное в некотором интервале около t 0 Если, кроме того, в этом интервале выполнено условие Липшица || f(t,x 1 ) – f(t,x 2 ) || L|| x 1 – x 2 || причем L не зависит от t, x 1 и x 2 , то это решение единственно и является непрерывной функцией от x 0 Для оценки значения постоянной Липшица L разложим в ряд Тэйлора функцию f(t,x 2 ) относительно точки x 1 (x 1 , x 2 R n ). 1 2 1 2 1 1 2 ( , ) ( , ) ( , )( ); ; [ , ]. i i i n f f t x f t x t x x x x x 106 Здесь 1 1 1 1 ( , ) n n n n f f x x f t x f f x x матрица Якоби. Таким образом, 2 1 2 1 || ( , ) ( , ) || || ( , ) || *|| ( ) || . f f t x f t x t x x x Все численные методы связаны с переходом от дифференциальной системы к разностной схеме и вычислением приближенных решений в форме сеточной функции. В дальнейшем будем обозначать x(t k ) решение дифференциальной системы в узле t k , а x k решение разностной схемы в том же узле, t k = t 0 + kh, где h шаг сетки, k = 0, 1, 2, … . Рассмотренный ранее аналитический метод решения задачи Коши применим только в случае линейных обыкновенных дифференциальных уравнений. 10.2 .Метод аналитического продолжения ( метод рядов). Предполагаем, что x, f(t,x) R n Предположим, что установлен факт корректности задачи и функция f(t,x) непрерывно дифференцируема требуемое число раз. Тогда решение x(t) можно разложить в ряд Тэйлора около центральной точки t k 107 2 1 ( ) ( 1) 1 ( 1) ( 1) ( 1) 1 ( ) ( ) '( ) ''( ) ... ( ) ( ) 2! ! ( 1)! ( ) [ ( , ( ), ) ( , ( ), )], ( , ( ), ) ( ), [ , ]; ( 1)! ( , ( ), ) '( ) ''( 2! s s s s k k k k k k s k k k k k k s s s k k k k k k k k k k k h h h x t x t hx t x t x t x s s x t h t x t h r t x h h r t x h x t t s h t x t h x t x t 1 ( ) ( 1) 1 ) ... ( ); ! , эквивалентное разностное уравнение имеет вид ( ) ( ) ( , ( ), ) ( ( ), ) (*) s s k s k k k k k k h x t s Следовательно x t x t t x t h r x h h ( 1) ( 1) ( ( ), ) ( ). ( 1)! s s s k k k h r x h x s Разностное уравнение метода (формулу метода) мы получим, если пренебрежем погрешностью r k 1 ( , , ) (**) k k k k x x t x h h 1 ( ) ( , , ) ' '' 2! ! s s k k k k k h h t x h x x x s Применение метода связано с использованием функции φ, определяемой на каждом шаге через производные до s – ого порядка включительно. Это предполагает использование соотношений: для производной первого порядка x'(t k )=f(t k ,x k ); для производной второго порядка x''(t k )= , ( ) , ( ) '( ) k k t t x x t k k k t t x x t f f x t t x и т. д. с последующим усложнением выражений. Здесь под , ( ) t t x x t k k f x понимается матрица размером n ×n, называемая матрицей Якоби 108 1 1 1 1 n n n n f f x x A f f x x Остается открытым вопрос о правомерности пренебрежения остаточным членом r k Назовем глобальной погрешностью величину E(h)= max|| x( t k )- x k || (k=1,2,..). Численный метод решения задачи Коши называют сходящимся, если для него E(h) → 0 при h→ 0. Принято говорить, что метод сходится с p-ым порядком точности (или имеет p-ый порядок точности), если для погрешности справедлива оценка E(h) Ch p Очевидно, чтобы быть работоспособным, метод должен быть сходящимся. Малая сеточная функция r k , добавленная к правой части уравнения (**), преобразует его в разностную систему (*). Её решение совпадает в узлах t k с решением исходного дифференциального уравнения. Величина r k выполняет роль погрешности аппроксимации. Говорят, что разностное уравнение метода (**) аппроксимирует исходное дифференциальное уравнение, если при h →0 r k →0. В нашем случае в результате выполнения этой операции уравнение (*), а значит и уравнение (**) совпадут с исходным дифференциальным уравнением, т. е. аппроксимация имеет место. Погрешность аппроксимации ( 1) ( 1) ( ( ), ) ( ). ( 1)! s s s k k k h r x h x s пропорциональна h s . В таком случае говорят, что уравнение метода (**) аппроксимирует дифференциальное уравнение с порядком точности s. 109 Для того чтобы разностная схема метода была сходящейся, остается доказать её устойчивость, т. е. ограниченную чувствительность к изменениям исходных параметров (в данном случае правой части). Действительно, пусть имеет место факт аппроксимации с порядком точности s : max k || r k || C 1 h s , ( 1) 1 1 ( ). ( 1)! s k C x s Добавляя к правой части разностной системы (**) r k , получим возмущенную систему(*), решение которой совпадает с решением исходного дифференциального уравнения (x(t k ), k=0, 1, 2, …) Если разностная система устойчива, то ||x(t k ) – x k || K max k || r k || KC 1 h s = C h s , k=1,2,… То есть max || ( ) || , C не зависит от k и h. s k k k E x t x Ch Иногда удобнее использовать л о к а л ь н у ю погрешность L. Локальной погрешностью называется погрешность на очередном k- ом шаге в предположении, что значение x k рассчитывается по формуле метода подстановкой в неё точных решений на предыдущих шагах. Из формул (*) и (**) следует L k+1 (t k ,x (s+1) (ξ k )) = h r k = α k h s+1 (в (**) заменяем x k на x(t k ) и затем вычитаем из(*); в результате получаем L k+1 =x(t k+1 )-x(t k )=hr k =α k h s+1 ). 10.2.1. Методы Эйлера. На практике находят применение следующие две модификации метода первого порядка ( в разложении учитываются слагаемые до первой производной включительно). Явный метод Эйлера x k+1 = x k + hf(t k ,x k ). Неявный метод Эйлера x k = x k+1 - hf(t k+1 ,x k+1 ). x k+1 = x k + hf(t k+1 ,x k+1 ). 110 Для определения x(t k ) в последней формуле нужно использовать разложение x(t) около центральной точки t k+1 Оценим устойчивость решений, получаемых этими методами. Как и прежде возмущенная система получается в результате изменения начальных условий. 0 0 * 0 0 . система ( ( ), ), нач. усл. x(t ) * . система ( * ( ), ), нач. усл. x*(t ) dx исх f x t t x dt dx возм f x t t x dt Обозначим v =x – x*. Вычитая из исходной системы возмущенную, получим 0 0 0 ( ( ), ) ( * ( ), ), нач. усл. v * dv f x t t f x t t x x dt Воспользуемся формулой разложения функции нескольких переменных в ряд Тэйлора с остаточным членом: f(x(t),t) - f(x*(t),t) = = * (( ), ) f x v t x (x(t) – x*(t)). Здесь * (( ), ) f x v t x матрица Якоби; ϴ R n×n , ϴ=diag( ϴ i ), 0 ϴ i 1, i= 1, 2 ,…, n. Тогда * 0 0 0 (( ), ) ( ); v * dv f x v t v t x x dt x . После подстановки в * (( ), ) f x v t x = A R n×n начальных условий это уравнение становится однородным дифференциальным уравнением с постоянными коэффициентами. Его используют в качестве тестового при анализе устойчивости. Таким образом, система тестовых уравнений в этом случае имеет вид v'(t) = Av(t). Пусть исходная задача асимптотически устойчива, то есть Re(λ i (A)) < 0; i=1,2,…,n. Естественно потребовать, чтобы система разностных уравнений метода также была асимптотически устойчивой. 111 Для явного метода Эйлера v k+1 =v k +hAv k =(E+hA)v k ; собственные числа λ (E+hA) =1+hλ(A). Условия асимптотической устойчивости |1+hλ i (A)| < 1, i=1,2,…,n. Иными словами [1+hRe λ i ] 2 +[hIm(λ i )] 2 < 1. Область устойчивости лежит внутри круга (см. рисунок). Таким образом, возникают ограничения на величину шага h < -2 Re( λ i )/( Re( λ i )^2 + Im(λ i )^2) . Так как тестовая задача предполагалась ассимптотически устойчивой (Re( λ i ) <0), то h < 2 |Re( λ i )|/( Re( λ i )^2 + Im(λ i )^2) . Если λ i вещественно, h < 2 |Re( λ i )|/ Re( λ i )^2 =2 /|Re( λ i )|. Если исходная задача неустойчива (Re(λ i (А))> 0 ), условия асимптотической устойчивости метода перестают выполняться при любом h > 0. 112 Решая тестовую задачу с использованием неявного метода Эйлера, получаем v k+1 =v k +hAv k+1 =(E-hA) -1 v k . Условия асимптотической устойчивости метода выполнены, если |1- λ i (A) | -1 < 1, т. е. |1- λ i (A) | > 1. Иными словами [1- hRe λ i ] 2 +[hIm(λ i )] 2 >1. Область устойчивости лежит вне круга (см. рисунок). Ограничение на шаг описывается формулой h >2 Re( λ i )/( Re( λ i )^2 + Im(λ i )^2) . В случае асимптотической устойчивости дифференциальной системы ( Re( λ i ) < 0 ) метод асимптотически устойчив при любом h. |