Диссертация (1145311), страница 30
Текст из файла (страница 30)
Учитывая, что наибольшая относительная погрешность вычисления exp(α) равна exp(∆(α)), где ∆(α) обозначает абсолютнуюпогрешность α, получаем следующую относительную погрешность вычисленияx∞ (V ):δ(x∞ (V )) = 3ε/2 + exp(|(V − θx )/dx |) − 1.217Рис. 4.6. Зависимость V от времени.Следовательно, искомая относительная погрешность равнаη = 12ε +7ε5ε++ 4 (exp(|ε(V − θm )/dm |)) +2 − 7ε/2 2 − 5ε/2+ exp(|ε(V − θn )/dn |) + exp(|ε(V − θh )/dh |).Это верхняя граница погрешности.На рис. 4.4 приведен потенциал V , найденный двумя различными методами (a: метод Эйлера с постоянным шагом интегрирования 0.1 мс; b: метод,приведенный в диссертации).
Программы для расчета с помощью обоих методовнаписаны на языке C++, вычисления проводились в арифметике с плавающейточкой с одинарной точностью. Время вычислений равно 0.047 с и 0.765 с соответственно. Относительная погрешность для одного шага метода вычисленияV не превосходит 0.25 и 0.0000255 соответственно.218Для тех же начальных данных (исключая V (0), которое было взято равным −50 мВ) та же система была решена с периодическими импульсными воздействиями тока (рис. 4.5). График зависимости напряжения от времени дляэтого случая представлен на рис. 4.6.
На рис. 4.7 показано, как влияет изменение Istim на синхронизирующие переменные.Рис. 4.7. Зависимость синхронизирующих переменных m, n, h от времени.При нахождении m максимальный шаг по времени не более 1.25·10−5 , наибольшая погрешность интегрирования на одном шаге не превосходит 10−5 . Принахождении V максимальный шаг по времени не больше 1.1 · 10−4 , наибольшаяпогрешность интегрирования на одном шаге не превосходит 10−5 .
Время счетасоставило 0.714286 с.Наибольший шаг равен 1.25419 мс, наименьший — 0.00158842 мс. Числошагов для интервала времени в 1000 мс оказалось равным 35719.219Шаг интегрирования вычислялся в каждой следующей точке. Следовательно, произведенные вычисления были сделаны с максимально возможнойточностью.Замечание 48. Все вычисления были произведены на компьютере со следующими спецификациями: AMD Sempron 3400+1.81 GHz. Для примеров времявычислений представлено для Virtual MS-DOS mode, 16-bit-application.
Последняя система уравнений была проинтегрирована для 32-bit Win-app.4.5. Оптимальное число шагов метода ЭйлераЗдесь мы рассмотрим модификацию предложенного алгоритма, котораяоснована на рассмотрении линейных систем обыкновенных дифференциальныхуравнений с постоянными коэффициентами.Рассмотрим задачу Коши для системы линейных дифференциальных уравненийẊ = AX,X(t0 ) = X0 (X0 6= 0),(4.27)здесь A = [aij ]mi,j=1 — вещественная квадратная матрица порядка m с постоянными элементами,X(t) = x01x1 (t)..
, X = .. .. 0xm (t)x0m.Будем сначала считать, что все элементы матрицы A и вектора X0 заданыточно.Решение задачи (4.27) имеет вид X = eA(t−t0 ) X0 . Найдем полную относительную погрешность приближенного решения задачи Коши методом Эйлера вточке t1 = t0 + 1.При данном числе шагов n, применяя метод Эйлера в точной арифметике,220получим следующее приближенное значение решения:n1X̃(t1 ) = E + A X0 .n(4.28)Обозначим относительную погрешность вычисления j-той компоненты вектора X̂(t1 ) по формуле (4.28) через δ1 (n).Норма вектора полной относительной погрешности, компоненты которогоравны полным относительным погрешностям компонент вектора X̂(t1 ), находится по формуле:m X xj (t1 ) − x̃j (t1 ) + εmn.F (n)=x̃j (t1 )(4.29)j=1Здесь первое слагаемое представляет собой норму вектора, состоящего из относительных погрешностей компонент найденного методом Эйлера решения вточной арифметике.
Второе слагаемое — норма вектора относительных вычислительных погрешностей соответствующих компонент решения. Таким образом, требуется найти значение n, при котором достигается минимум функцииF (n).Рассмотрим общий случай: ни одна из компонент вектора X̂(t1 ) не обращается в нуль.Определение 40. Назовем оптимальным числом шагов метода Эйлера число шагов nopt , при котором значение решения в точке находится с наименьшей полной относительной погрешностью.Оптимальное число шагов метода Эйлера можно найти, воспользовавшисьследующей теоремой.Теорема 71.
Оптимальное число шагов nopt для нахождения значения решения системы дифференциальных уравнений (4.27) при t1 = t0 +1 приближенноможет быть найдено по формулеv u ẍm (t1 ) u ẍ1 (t1 ) ẍ2 (t1 ) ut x1 (t1 ) + x2 (t1 ) + · · · + xm (t1 ) nopt ≈.2mε(4.30)221Доказательство. Пусть AJ — жорданова нормальная форма матрицыA, а T = [tij ]mi,j=1 — матрица, составленная из векторов канонического базиса,такая, что AJ = T −1 AT .Тогда равенство (4.28) можно записать в видеn1X̂(t1 ) ≈ T E + AJ T −1 X0 .n(4.31)Рассмотрим сначала случай различных собственных чисел матрицы A.Пусть A — квадратная матрица порядка m.
Предположим, что матрицаA имеет l пар комплексно-сопряженных собственных чисел µ1 , µ̄1 , . . . , µl , µ̄l иm − 2l вещественных собственных числа λ1 , λ2 , . . . , λm−2l . (Эти собственныечисла считаем точными, однако их нахождения не требуется.) Для комплексных собственных чисел будем также использовать тригонометрическую формуµj = ρj (cos ζj + i sin ζj ), j = 1, 2, . . . , l, причем будем считать sin ζj > 0.Тогда согласно (4.31) j-я компонента вектора X(t1 ) — решения задачиКоши (4.27) будет иметь вид(j)(j)(j)(j)(j)(j)(j)c1 eλ1 + c2 eλ2 + · · · + cm−2l eλm−2l + d1 eµ1 + d¯1 eµ̄1 + · · · + dl eµl + d¯l eµ̄l ,(j) (j) (j)где постоянные коэффициенты ck , dk , d¯1 зависят от элементов матрицы T .Соответственно, j-я компонента вектора X̃(t1 ) — приближенного решениязадачи Коши методом Эйлера имеет видnnλ1λm−2l(j)(j)c1 1 ++ · · · + cm−2l 1 ++nnµ1 n ¯(j) µ̄1 nµl n ¯(j) µ̄l n(j)(j)+d1 1 ++ d1 1 ++ · · · + dl 1 ++ dl 1 +.(4.32)nnnnm (j) λ1(j)(j) µ1(j) µ̄1(j) µl(j) µ̄lλ¯¯m−2lX c1 e + · · · + cm−2l e+ d1 e + d1 e + · · · + dl e + dl e + εmn−1F (n) =P Pm−2l (j)nn(j)(j)lµk nµ̄λkk1 + n + k=1 dk 1 + n + d¯k 1 + nj=1 k=1 ckилиm−2llX (j)X(j) ρk cos ζk(j)λkcosρsinζ+χce+2rekkkkkm Xk=1k=1F (n)=−1 m−2l +εmn,n Xn/2lX2λρρcosζj=1 kkk(j)(j)(j)k+2rk 1+ 2 +2cos ωk n+χkck 1+nnn k=1k=1222где(j)dk=(j)rk(j)cos χk+(j)i sin χkζk1 + ρk cosµk n= arccos q,, ωk = arg 1 +nρ2kρk cos ζk1 + n2 + 2 nс учетом положительности sin ζk , k = 1, 2, .
. . , l.Для доказательства теоремы нам понадобятся следующие формулы:n 2λk1(j)(j) λk(j) λk λkck 1 ++o,= ck e − ck en2nnµk n ¯(j) µ̄k n(j)(j) ρk cos ζk(j)dk 1 ++ dk 1 += 2rk ecos ρk sin ζk + χk −(4.33)nn 1 (j) 2 ρk cos ζk1(j)− r k ρk ecos (ρk sin ζk + 2ζk + χk + o.nnРазложим F (n) в сходящийся ряд по степеням 1/n до членов второго порядка, используя равенства (4.33), получим функцию F̂ (n): m−2llXX(j) 2 λk(j) 2 ρk cos ζk(j) c k λk e +2rk ρk ecos ρk sin ζk + 2ζk + χk m X1 k=1k=1F̂ (n) =+εmn.lm−2l2n j=1 XX(j)(j) ρk cos ζk(j) λkecosρsinζ+χe+2rckkkkkk=1k=1Продифференцируем F̂ (n) по n и приравняем полученное выражение к нулю,получим приближенное уравнение для нахождения n:vPuPl(j)(j)(j)m−2l22λρcosζmuXk cosρk sin ζk + 2ζk + χk k=1 ck λk e k + k=1 2rk ρk e ku /2mε.n≈tPm−2l (j) λPl(j)(j)ρk cos ζk cos ρ sin ζ + χk +ce2rej=1 kkk=1k=1kkk(4.34)Осталось заметить, чтоPm−2l (j) 2 λk Pl(j) 2 ρk cos ζk(j)cos ρk sin ζk + 2ζk + χkk=1 ck λk e +k=1 2rk ρk eẍj (t) =.Pm−2l (j) λPl(j) ρk cos ζk(j)xj (t) t=t0 +1kc e +2r ecos ρ sin ζ + χk=1kk=1kТем самым, справедливость формулы (4.30) установлена.kkk223Теперь предположим, что у матрицы A имеются кратные собственные числа.
В этом случае нам понадобятся формулыn−pk (n−1)(n−2) . . . (n−pk +1)λkλ2k +2pk λk +pk (pk −1)eλk11+1−+o=pk !npk −1npk !2nn(4.35)и(j) (n−1)(n−2) . . . (n−pk +1)dkpk !npk −1µk n−pk1++n µ̄ n−pkk(j) (n−1)(n−2) . . . (n−pk +1)¯1+=+ dkp−1pk !n kn(j)(4.36)(j)rk ρk cos ζk h 22rk ρk cos ζk(j)(j)=ecos(ρk sin ζk + χk ) −eρk cos(ρk sin ζk + χk + 2ζk ) +pk !npk ! i1(j)(j)+ 2pk ρk cos(ρk sin ζk + χk + ζk ) + pk (pk − 1) cos(ρk sin ζk + χk ) + onв обозначениях, введенных ранее. Аналогично рассмотренному ранее случаю,получаем, что формула (4.30) остается справедливой и в случае наличия кратных собственных чисел у матрицы A.Осталось проверить, что производныефункций F (n) и F̂ (n) достаточно 1близки.
Для этого рассмотрим o0для произвольного простого собственногоnчисла ν матрицы A: 3 0 1 8ν + 3ν 41 νo = |ce | +o. 12n3n n4 1 =Для кратного собственного числа ν кратности p имеем o0n 432232ν −3ν −4(3p+2)ν −3p(5p+3)ν +6p(1−p )ν−p(3p −10p +9p−2)+= |ce | 12p!n3 1 +o,n4 откуда сразу же следует близость производных функций F (n) и F̂ (n) при больших n.224 1Одновременно мы оценили скорость убывания o. Она пропорциональnна 1/n3 .2Замечание 49.
Заметим, что нормальная форма Жордана матрицы A невычисляется и нигде в алгоритме она не используется. Она применяетсятолько для сравнения точного результата, который находится с ее помощью, с результатом, найденным приближенно.Пусть sj (t) = sign ẍj (t)/xj (t) , j = 1, 2, . . . , m.Так как Ẍ(t) = A2 X(t), справедливо следующее утверждение.Следствие 22. Оптимальное число шагов nopt для нахождения значения решения системы дифференциальных уравнений (4.27) при t1 = t0 + 1 приближенно может быть найдено по формулеvuux (t )u 1 1u u x2 (t1 )sm (t1 )s1 (t1 ) s2 (t1 )u 1nopt ≈ u,,...,A2 u 2mε x1 (t1 ) x2 (t1 )xm (t1 ) ...utxm (t1 ).Теперь предположим, что элементы матрицы A заданы с относительнымипогрешностями, не превосходящими δA, а элементы вектора X0 — с относительными погрешностями, не превосходящими δX0 .В этом случае относительная вычислительная погрешность вычисленийдля любой компоненты вектора X̃(t1 ) равна (ε + δA)n + δX0 .Следствие 23.
Оптимальное число шагов nopt для нахождения значения решения системы дифференциальных уравнений (4.27) при t1 = t0 + 1 прибли-225женно может быть найдено по формулеvuuu x1 (t1 )u s1 (t1 ) s2 (t1 )sm (t1 )2nopt ≈ (1 + δA)u...u x (t ) , x (t ) , . . . , x (t ) A 2 1m 1t 1 1xm (t1 ) /2m(ε + δA),(4.37)в обозначениях следствия 22, где δA — максимальная относительная погрешность элементов матрицы A.Доказательство следствия 23 аналогично доказательству теоремы 71.Теорема 72. Для оптимального числа шагов nopt относительная погрешностьметода Эйлера на шаге совпадает с относительной вычислительной погрешностью на шаге.Доказательство.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
