Диссертация (1145311), страница 28
Текст из файла (страница 28)
е. существует константа L ≥ 0 такая, что для всех Y 6= Xвыполнено неравенство||F (t, X) − F (t, Y )|| ≤ L||X − Y ||.(4.9)200Предлагаемый метод позволяет вычислить решение системы дифференциальных уравнений достаточно точно. Метод использует интегрирование с оптимальным шагом. Как это видно из приведенных далее численных примеров,данный метод может быть применен в том числе и для жестких систем. Времявычислений довольно мало, что очень важно для расчетов в режиме реальноговремени. Один из примеров работы предложенного метода представляет собойрасчеты для модели ионных каналов — в этом случае как раз применяютсявычисления в режиме реального времени, так как компьютер должен обрабатывать результаты эксперимента в процессе самого эксперимента.Далее используются соответствующие векторная и матричная нормы:x 1 Xm×1 = ...
: ||X|| = |x1 | + |x2 | + · · · + |xm |,xm||B|| = max ||Bj ||,1≤j≤mгде B = (B1 , B2 , . . . , Bm ) — квадратная матрица размерности m × m со столбцами B1 , B2 , . . ., Bm .Известно, что в конечномерном пространстве все нормы эквивалентны.Следовательно, выбор нормы может быть произведен из соображений удобства.В дальнейшем через X(t) будем обозначать точное решение задачи Коши (4.8), X̃(t) — решение, полученное с помощью метода Эйлера и учитывающее только ошибку метода, и X̂(t) — вычисленное решение. Будем рассматривать полную относительную погрешность приближенного решения задачи Коши в точке t0 + h.
Эта погрешность является суммой погрешности метода иошибок округления.В методе Эйлера есть два источника ошибок округления. Первая ошибка возникает при вычислении F (t, X). Норму этой ошибки обозначим через η.Верхняя граница данной ошибки может быть получена с помощью примененияпринципов вычисления в арифметике с плавающей точкой [61, 99]. Вторая ошиб-201ка получается из-за ошибок округления при применении формулы Эйлера.
Ееверхнюю границу можно найти, используя результат из книг, приведенных выше: = εm.Как уже говорилось в разделе 2.4.2 главы 2, точность вычислений можетбыть увеличена, если использовать результат, представленный в книге [61].Вычисления с большей точностью, чем двойная точность, выполняются спомощью специальных программ, т. е. они зависят от компьютера. Такие программы трудно сделать быстрыми и машинонезависимыми. Однако, достаточнохорошие результаты можно получить, производя вычисления с одинарной точностью (float).Для одного шага метода Эйлера h очевидно имеемX̃(t0 + h)= X0 + F (t0 , X0 )h,X̂(t0 + h)= X0 + F (t0 , X0 )h ± (εm + η)X(t0 + h) == (X0 + F (t0 , X0 )h)(1 ± (εm + η)).Определение 38. [93] Выражение 1 (X(t + τ ) − X(t)) − F (t, X(t)) τназывается локальной ошибкой дискретизации метода Эйлера в точке t.Эта ошибка является мерой того, насколько отличается Ẋ(t) от F (t, X(t)).Определение 39.
Ошибка дискретизации в точке t + τ равна величине||X̃(t0 + τ ) − X(t0 + τ )||.Пусть h1 , h2 , . . . , hn — n шагов метода Эйлера. Пусть Hk обозначает суммуh1 + h2 + . . . + hk , Hn = T .Максимальное значение ошибки дискретизации на интервале [t0 , t0 + T ],т. е.E(t0 , T ) = max ||X̃(t0 + Hi ) − X(t0 + Hi )||,1≤i≤N202называется глобальной ошибкой дискретизации, или глобальной погрешностьюдискретизации на промежутке [t0 , t0 + T ].Верхняя граница нормы абсолютной погрешности может быть выраженачерез глобальную ошибку дискретизации (см., например, [93]):eLT − 1||X(t + τ ) − X̃(t + τ )|| ≤M τ,2L(4.10)если τ достаточно мало (iτ ≤ T ) иX(t+τ)−X(t) ≤ M τ.F (t, X) −τТеперь рассмотрим полную погрешность, которая равна сумме ошибки дискретизации и ошибок округления.
Далее мы будем считать, что выполняетсяпредположение nυ < 0.1, которое на практике действительно справедливо вбольшинстве случаев. Следующая теорема показывает, что оптимальный шагинтегрирования дает минимальное значение полной погрешности, т. е. он является в самом деле оптимальным в смысле вычислительной точности.Теорема 69. Полная погрешность минимальна, когда погрешность методасовпадает с ошибками округления.Доказательство. Рассмотрим решение задачи Коши в точке t1 = t0 + T .Предположим, что h1 , h2 , . . . , hn достаточно малы.
Через h обозначим max hk .k=1,...,nТогда норма абсолютной ошибки дискретизации решения в точке t1 такова (4.10),что выполняется следующее условие:||X(t1 ) − X̃(t1 )|| ≤Положим Fmax =eLT − 1M h.2Lmax ||F (t, X(t))|| и оценим ошибку округления решеt∈[t0 ,t0 +T ]ния в арифметике с плавающей точкой.Выполняется следующее равенство:X̂(t0 + Hk ) = (X̂(t0 + Hk−1 ) + hk F (t0 + Hk−1 , X̂(t0 + Hk−1 )(1 ± η))(1 + δk ),203где 0 ≤ |δk | ≤ υ. Для вычисленного решения в точке t1 имеем||X̂(t1 ) − X̃(t1 )|| ≤ ||X0 (1 + δ1 ) . .
. (1 + δn )++F (t0 , X0 )h1 (1 ± η)(1 + δ1 ) . . . (1 + δn )++F (t0 + H1 , X̂(t0 + H1 ))(1 ± η)h2 (1 + δ2 ) . . . (1 + δn ) + . . . ++F (t0 + Hn−1 , X̂(t0 + Hn−1 ))(1 ± η)hn (1 + δn ) − X̃(t1 )||.Следовательно, верхняя граница нормы абсолютной погрешности для X̂(t0 +Hn ) удовлетворяет неравенству:γn (||X0 || + T Fmax (1 + η)) ≤ 1.06nυ (||X0 || + T Fmax (1 + η)) ,nυ(см. [99], где доказана последняя оценка).1 − nυ/2Из неравенства между средним арифметическим и средним геометриче-где γn =ским получаем, что верхняя граница полной погрешности минимальна, когдавыполняется следующее равенство:eLT − 1M h = 1.06nυ (||X0 || + T Fmax (1 + η)) ,2L(4.11)т. е.
верхняя граница полной погрешности минимальна, когда верхние границыошибки дискретизации и ошибки округления совпадают.Приведенный результат выполняется для произвольного числа шагов метода Эйлера. Следовательно, он справедлив и для одного шага (разумеется,шаги метода должны быть достаточно малыми).
Тем самым, теорема доказана.Теорема 69 может быть использована при минимизации полной погрешности. Если при численном интегрировании системы дифференциальных уравнений каждый шаг метода Эйлера будет оптимальным, то полная погрешностьрешения будет наименьшей.204К сожалению, не существует способа найти оптимальный шаг метода Эйлера для уравнения (4.11) без дополнительных вычислений.
Для нахожденияего воспользуемся следующей теоремой:Теорема 70. Локально оптимальный шаг метода Эйлера может быть найден из следующего уравненияhoptvu2(η + εm)u .= u t ẍj (t0 +hopt ) m xj (t0 +hopt ) j=1 (4.12)Доказательство. Для одного шага метода Эйлера покомпонентная ошибка дискретизации решения задачи Коши удовлетворяет соотношениюdh2X(t0 + h) − X̃(t0 + h) = F (t0 + θ, X(t0 + θ)) ,dt2(4.13)где 0 < θ < h.Следовательно, относительная ошибка дискретизации метода для одногошага равна ẍ (t + θ) m h2 ẍ (t + h) m h2 j 0 j 0||δX(t0 + h)|| = ≈ , xj (t0 + h) j=1 2 xj (t0 + h) j=1 2если шаг метода h достаточно мал.В самом деле, имеемd3 XẌ(t0 + h) = Ẍ(t0 + θ) + (h − θ) 3 (t0 + Θ),dtгде Θ принадлежит отрезку с концами в точках 2 t0 + h и t0+ θ. Таким образом, d Fесли h меньше, чем ε/M, где M = max 2 (t, X(t)), нормы Ẍ(t0 + h) иdtt∈[t0 ;t0 +h]Ẍ(t0 + θ) отличаются меньше, чем на ε.Эта погрешность для оптимального шага интегрирования должна совпадать с ошибкой округления, откуда получаем ẍ (t + h ) m h2 opt j 0opt= η + εm.
xj (t0 + hopt ) j=1 2205Следовательно, оптимальный шаг метода Эйлера равенvu2(η + εm)u ,hopt = u t ẍj (t0 +hopt ) m xj (t0 +hopt ) j=1 что и завершает доказательство.Следствие 21. Если начальные данные задачи Коши имеют погрешность,то справедлива следующая формула:vu 2(ε + εm + η)u0 .hopt = u t ẍj (t0 +hopt ) m xj (t0 +hopt ) j=1 Здесь ε0 обозначает относительную погрешность начальных данных.4.3. Реализация методаРассмотрим простой итеративный метод для нахождения локально оптимального шага метода Эйлера. Сначала требуется найти приближенное значение X(t0 + h1 ) для достаточно малого шага h1 .
Затем находим новый шаг иновое приближенное значение решения с помощью следующих формул:vu 2(εm + η)u , X̃ = X0 + hk F (t0 , X0 ),hk+1 = u (4.14)t ẍj (t0 +hk ) m k x̃j (t0 +hk ) j=1 где k = 1.Известно [93], что ошибка дискретизации в методе Эйлера стремится к нулю, если шаг метода также стремится к нулю. Из формул (4.14) следует, чтоесли шаг hk меньше оптимального, то шаг hk+1 будет больше (hk+1 > hk ) и наоборот (если hk больше оптимального шага, то hk+1 < hk ).
Предположим, что h1выбран достаточно близко к hopt . Тогда, учитывая то, что что все вектор-функции и функции в уравнении (4.14) непрерывны, получаем, что чем меньше hk206отличается от hopt , тем меньше становится норма относительной погрешностиX̃(t0 + hk ). Обратное тоже верно.Следовательно, нужно доказать, что возможно выбрать hk достаточноблизко к hopt .Сначала положим шаг h1 равным ε (или другому минимально возможному значению). Это нижняя граница шага интегрирования.
Используя формулы (4.14), вычислим новый шаг интегрирования. Если его величина получитсяменьше ε, положим его равным ε и найдем значение решения в новой точкеt0 + h1 .Если новый шаг интегрирования h2 будет больше ε, вычислим шаг интегрирования h3 по формуле (4.14). Если h3 < h2 , шаг h2 является верхней границейдля оптимального шага интегрирования. В противном случае будем вычислятьвеличины h4 , h5 и т.д., до тех пор, пока не получим hj+1 ≤ hj . Тогда для некоторой итерации получим hj = hj+1 = hopt (в том случае, если последовательностьhj возрастает, и ее верхняя граница равна hopt ) или hj+1 < hj .Таким образом, мы получаем текущий шаг интегрирования hk и границыоптимального шага интегрирования (hmin и hmax ).