Диссертация (1145311), страница 31

Просмтор этого файла доступен только зарегистрированным пользователям. Но у нас супер быстрая регистрация: достаточно только электронной почты!

Текст из файла (страница 31)

Воспользуемся формулами (2.21) для случая различных собственных чисел матрицы A и формулами (4.35), (4.36) соответственнодля кратного вещественного собственного числа λk кратности pk и для кратныхкомплексных собственных чисел µk и µ̄k кратностей pk (для случая недиагональной формы Жордана).Пусть δξ — вектор относительных погрешностей решения на одном шаге,то есть j-я компонента δξ равна отношению j-й компоненты вектора X̄(t1 ) −eA/n X0 к j-й компоненте X̄(t1 ).С учетом (2.21), (4.35), (4.36) приближенное равенство||δξ|| ≈ mεравносильно приближенному равенству (4.30).2Следствие 24. При нахождении значения решения задачи Коши (4.27) в точке t1 = t0 + τ методом Эйлера имеемnAτX̃(t1 ) = E +X0 .n226Таким образом, этот случай сводится к случаю τ = 1 заменой матрицы A наматрицу Aτ .

Следовательно, достаточно рассмотреть только случай τ = 1.Из формулы (4.37) и замечания 24 следует, что погрешность элементовматрицы A влияет не только на значение решения в точке, но и на саму точку,в которой вычисляется решение. Действительно, поскольку матрица системыстановится равной A(1 ± δA), то вместо точки t0 + 1 значение решения вычисляется в точке t0 + 1 ± δA.4.5.1. Вычислительный алгоритмПредложенный здесь алгоритм позволяет не только найти оптимальноечисло шагов метода Эйлера, но и вычислить значение решения задачи Кошимаксимально точно.Воспользуемся методом Банаха (простых итераций).

Сначала найдем приближенное значение X1 (t1 ), которое получается интегрированием с помощьюметода Эйлера с достаточно большим числом шагов n1 . Затем по найденномузначению X1 (t1 ) определим следующее значение n2 числа шагов метода Эйлера, и построим новое решение с числом шагов n2 . И т.д. Каждое следующеезначение nk+1 числа шагов будем находить по формулеnk+1v!u (k)(k)u s (t1 )s(t)m1t 12 X (t )/2mε , X (t ) = =,...,Ak 1k 1(k)(k)x1 (t1 )xm (t1 )(k)x1 (t1 )(k)x2 (t1 )...(k)xm (t1 ).(4.38)Здесь для вещественного α величина dαe обозначает наименьшее целоечисло, которое больше или равно α.Теорема 73.

Последовательность (4.38) сходится, если n1 достаточно велико.227Доказательство. Обозначим через ∆n = nk − nopt , и черезT∆xm∆x1δXpoln =,...,,x1xmгде ∆xj = xj (nk ) − xj (nopt ) = xj (nk ) − xj . Также через A2 (j) и A2k (j, k =1, 2, . . . , m) обозначим соответственно j-ю строку и k-й столбец матрицы A2 .Далее будем пользоваться равенством1∆xj∆xj1=− 2 +o,j = 1, . . . , m.xj + ∆xjxjxjxj√√∆a∆aПоскольку a + ∆a − a = √ + o √ , то для оценки близости неко2 aaторого n и nopt рассмотрим разность  x1 + ∆x1   x1 s1sms1sm..,...,A2 ,...,A2  ...  =−.x1 + ∆x1xm + ∆xmx1xmxm + ∆xmxm2 (1)s /x A  x1  1 1 . ∆x1∆xm   ..  +..−,...,.x1xm xmsm /xm A2 (m)∆x1 /x1sms1.22 .,...,(x1 A1 , .

. . , xm Am ) + + o(||δXpoln ||)..x1xm∆xm /xmТогда  x1 + ∆x1  ss1smsm1.2..,...,A ,...,A2 − x1 + ∆x1xm + ∆xmx1xmxm + ∆xm2 (1)s /x A 1 1  x 1    .    ..  ||δXpoln || +.≤ (1, . . . , 1) ..xm sm /xm A2 (m)x1...xm ≤228 ssm 122 + ,...,(x1 A1 , . .

. , xm Am )  x 1xm1...1    ||δXpoln || + o(||δXpoln ||) = = 4mεn2opt ||δXpoln || + o(||δXpoln ||),так как(1, . . . , 1) s1 /x1 A2 (1)...sm /xm A2 (m)x  1  ss  .. 1m,...,A2  .  =x1xmxmx1..  = 2mεn2. optxmиsms1,...,x1xm1   ssm122  .. (x1 A1 , . . . , xm Am )  .  =,...,A2 x1xm 1Отсюда|nk − nopt | ≤√x1..  = 2mεn2 .. optxm2mεnopt ||δXpoln || + o(||δXpoln ||).Для оценки ||δXpoln || воспользуемся равенствами (2.21), (4.35), (4.36).Имеем1111 + 2mε + o||δXpoln || ≤ 2mεn2opt −+o.nk−1 nopt nk−1noptСледовательно, при достаточно большом nk с уменьшением относительнойпогрешности решения ||δXpolnk ||, следующее значение nk+1 становится ближе кnopt , а при приближении nk к nopt относительная погрешность приближенногорешения уменьшается.Следствие 25. Так как для правой части системы уравнений (4.27) выполняется условие Липшица с L = ||A||, а кроме того при τ ≤ 1 τX(t+τ)−X(t)AX(t) − ≤ ||A||2 e||A|| ||X0 ||, 2τ229согласно [8], оценим абсолютную погрешность метода Эйлера при t1 = t0 + 1.Она для данного числа шагов n по норме не превосходит||A||e||A|| ||A||||∆(n)|| ≤e− 1 ||X0 ||.2nТем самым, n можно выбрать так, чтобы ∆(n) ||A||e||A|| ||A||≤||δξ|| min |xj (t1 )| ≤ e− 1 ||X0 || ≤ mε min |xj (t1 )|.n 2n2j=1,mj=1,mПосколькуmin |xj (t1 )| ≥ min |x0j |e−||A|| ,j=1,mj=1,mто если у вектора X0 нет нулевых компонент, при выборе n1 исходя из условияvu ||A|| e||A|| − 1 ||X ||u0n1 ≥ e||A|| t,2mε min |x0j |j=1,mмы заведомо получим n1 ≥ nopt .Однако последовательность (4.38) во многих случаях сходится и при выборе n1 = 1 (так, приведенный пример просчитан при данном значении n1 ).Для найденного оптимального значения числа шагов интегрирования приведем некоторую оценку погрешности вычисления значения решения задачиКоши методом Эйлера в точке.Теорема 74.

Для порядка матрицы A m < 42016 при точности float (ε =1.19·10−7 ) полная погрешность метода составит не более 1% при выполненииусловия(1)s1 (t1 ),...,(1)x1 (t1 )(1)sm (t1 )(1)xm (t1 )!A2 X1 (t1 ) < 2mε.Доказательство. Из теоремы 72 и условия утверждения следует, что приоптимальном значении nopt = 1 полная погрешность не превосходит 2mnε =0.01.230Замечание 50. В приведенном алгоритме рассмотрен общий случай, когдани одна из компонент вектора Xk (t1 ) не обращается в нуль. Если же в процессе &вычислений' какие-то компоненты обнуляются, то nopt положим равr||A2 ||ным.

Действительно, переходя в этом случае к абсолютной по2mεгрешности, применяя формулы (2.21), (4.35) и (4.36), аналогично результатутеоремы 71, получаемsnopt ≈||Ẍ(t1 )||.2mε||X(t1 )||Теперь воспользуемся свойствами нормы для получения верхней границы n:ssr||Ẍ(t1 )||||A2 X(t1 )||||A2 ||=≤.(4.39)2mε||X(t1 )||2mε||X(t1 )||2mεТем самым, при ненулевом векторе Xk (t1 ), который однако имеет по крайнеймере одну нулевую компоненту, мы можем положить'&r2||A ||nopt =.2mεНеравенством (4.39) можно воспользоваться для определения n1 в томслучае, когда вектора X0 имеет хотя бы одну нулевую компоненту. В этомслучае можно положить&rn1 =||A2 ||2mε'.Пример 15. В качестве еще одного примера приведем построение алгебраической кривой F (x, y) = 0, гдеF (x, y) = −3840x6 + 512x5 y − 2688x4 y 2 + 256x3 y 3 − 624x2 y 4 ++ 32xy 5 − 48y 6 + 7552x5 + 6784x4 y − 544x3 y 2 ++ 416x2 y 3 − 80xy 4 − 272y 5 + 1168x4 − 9792x3 y ++ 7280x2 y 2 + 288xy 3 − 1340y 4 − 7200x3 − 4896x2 y −− 216xy 2 − 2520y 3 − 120x2 + 5616xy − 4164y 2 ++ 2016x + 2016y + 441.231Будем строить данную кривую как решение системы обыкновенных дифференциальных уравнений.

Для этого рассмотрим следующую систему ОДУdx ∂F dy∂F=,=−.dt∂y dt∂xПерепишем эту систему в виде  xẋ  = A(x, y)   ,yẏ(4.40)где элементы матрицы A следующиеa11 (x, y) = 2016/x + 512x4 − 5376x3 y + 768x2 y 2 − 2496xy 3 + 160y 4 ++ 6784x3 − 1088x2 y + 1248xy 2 − 320y 3 − 9792x2 ++ 14560xy + 864y 2 − 4896x − 432y + 5616, если x 6= 0,a11 (x, y) = 160y 4 − 320y 3 + 864y 2 − 432y + 5616, если x = 0,a12 (x, y) = −288y 4 − 1360y 3 − 5360y 2 − 7560y − 8328, если x 6= 0,a12 (x, y) = 2016/y − 288y 4 − 1360y 3 − 5360y 2 − 7560y − 8328,если x = 0,a21 (x, y) = −2016/x − 2560x3 y + 10752x2 y 2 − 768xy 3 + 1248y 4 −− 27136x2 y + 1632xy 2 − 832y 3 + 29376xy − 14560y 2 ++ 9792y + 23040x4 − 4672x2 − 37760x3 ++ 21600x + 240, если x 6= 0,a21 (x, y) = 1248y 4 − 832y 3 − 14560y 2 + 9792y + 240, если x = 0,a22 (x, y) = −32y 4 + 80y 3 − 288y 2 + 216y − 5616, если x 6= 0,a22 (x, y) = −2016/y − 32y 4 + 80y 3 − 288y 2 + 216y − 5616,если x = 0.Эта система была проинтегрирована с помощью предложенного метода.Были использованы оба варианта применения рассмотренного алгоритма.232Первый вариант.

Фиксированы значения t, в которых матрица системыпересчитывается, а значение решения для данных t находится с помощью оптимального числа шагов.Второй вариант. В каждой точке t находится оптимальный шаг интегрирования и значение решения в новой полученной точке, если использовать этотоптимальный шаг. Матрица системы в каждой новой полученной точке пересчитывается.Чтобы не возникала ситуация overflow при написании программы на C++матрица системы была домножена на 10−8 , тем самым уменьшился исходныйшаг интегрирования. В связи с этим оптимальное число шагов практически вовсех случаях стало равно 1 и два предложенные варианта для системы уравнений (4.40) дали совпадающие решения, которые и представлены на рис. 4.8.Также на рис. 4.8 показана аналитическая кривая, которая выделена более жирно.

Как видно из этого рисунка, кривые отличаются незначительно.При интегрировании системы (4.40) использовалась таблица 2 начальныхданных задачи Коши.Замечание 51. Система дифференциальных уравненийdX= AXdt(4.41)заменой переменных t = kτ , k 6= 0 сводится к системеdX(kτ )dZ= kAX(kτ ) или= kAZ, (Z(τ ) = X(kτ )),dτdτ(4.42)т. е. системе такого же вида с матрицей kA.Найдя оптимальное число шагов метода Эйлера для системы (4.41) по формуле (4.14), мы найдем также пропорциональное ему оптимальное число шаговдля системы (4.42). Нужно только вычислить число шагов, которое получается, вообще говоря, не целым, с необходимой точностью. Поделив матрицу A нанайденное для нее оптимальное число шагов, получим матрицу системы, для233Рис.

Характеристики

Список файлов диссертации