Диссертация (1145311), страница 29
Текст из файла (страница 29)
Найдем hk+1 . Если hmin <hk+1 < hmax , мы можем указанным выше образом вычислить следующий шагинтегрирования и улучшить границы для оптимального шага.Если hk+1 < hmin или hk+1 > hmax , тогда положим hk+1 равным (hmin +hmax )/2. Очевидно что в этом случае новый шаг hk+1 станет либо верхней, либонижней границей для оптимального шага интегрирования.Мы останавливаемся, когда hk+1 = hk . (В арифметике с плавающей точкойэто эквивалентно выполнению условия |hk+1 − hk | < ε.)Каждая итерация дает все более точные границы для оптимального шага интегрирования.
Расстояние между hmin и hmax стремится к нулю. Такимобразом, за конечное число итераций мы получим оптимальный шаг методаЭйлера.207Замечание 38. Мы предполагали, что все компоненты вектора X̃k ненулевые. Если же хотя бы одна компонента данного вектора равна нулю, относительная погрешность не определена. В этом случае можно использоватьабсолютную погрешность. Тогда первая формула равенств (4.14) принимаетвидshk+1 =2(η + εm)||X(t0 + hk )||.||Ẍ(t0 + hk )||(4.15)С помощью (4.13) получаем следующее приближенное равенство для абсолютной погрешности: 2 2m hm h||∆X(t0 + h)|| = (ẍj (t0 + θ))j=1 ≈ (ẍj (t0 + h))j=1 ,22где h достаточно мало.Для нахождения оптимального шага метода Эйлера мы должны приравнять эту погрешность к норме абсолютной погрешности (ошибок округления) (η + εm)||(xj (t0 ) + h)mj=1 || для оптимального шага hopt .
Следовательно,получаем формулу (4.15).Замечание 39. Шаг метода Эйлера, приведенный в книге [106] был полученбез строгого доказательства, в некоторой степени интуитивно. Верхняя граница этого шага совпадает с величиной шага, которая может быть полученас помощью одной итерации по формуле (4.14), где h1 = 0.Действительно, для одного дифференциального уравнения ẋ = f (x, t) дляодного шага в книге [106] предлагается выбирать шаг h такой, чтоs!2ε|x(t0 + h)|h < min, hmax ,|ẍ0,1 |f (t0 + h, x(t0 + h)) − f (t0 , x0 ). Подставляя в это неравенство h1 =h0, m = 1, получаемs!2ε|x0 |h < min, hmax .|ẍ(t0 )|где ẍ0,1 =208f (t0 + h, x(t0 + h)) − f (t0 , x0 )(Здесь ẍ0,0 = lim= ẍ(t0 ).) Следовательно, граниh→0hs2ε|x0 |в правой части полученного неравенства в точности совпадает сца|ẍ(t0 )|границей, полученной по формулам (4.14).Сами формулы для нахождения оптимального шага интегрирования достаточно просты, однако в их применении возникает ряд сложностей.
Наиболеесложной является оценка погрешности η. Для того чтобы найти ее, можно воспользоваться подходом, представленным в книге [61].Замечание 40. Чтобы сократить время вычислений, производные, стоящиев правых частях системы, могут быть заменены конечными разностями, какэто предложено в книге [106].Очевидно, что нет необходимости менять шаг интегрирования при каждомпереходе к новой точке.Замечание 41. Ясно, что существует очевидная нижняя граница для шагаинтегрирования. Наименьший шаг интегрирования не может быть меньшеε. Причина этого в том, что εF (t, X(t)) дает абсолютную погрешность вычисления F (t, X).Замечание 42.
Как показано в Примере 3, метод может быть применени для системы дифференциальных уравнений с кусочно-дифференцируемымиправыми частями.Замечание 43. Во многих случаях (см., например, примеры 12 и 13), формулы (4.14) могут быть применены без оценок hmin and hmax , что довольносильно сокращает время работы процессора.Алгоритм. В каждой точке определяется оптимальный шаг интегрирования, исходя из условия, что ошибка округления совпадает с ошибкой дискретизации. Алгоритм следующий:2091.
Положим начальный шаг интегрирования h0 = ε, нижнюю и верхнююграницы для него положим равными hmin = ε, hmax = 1.0 (hmin и hmax могутбыть выбраны, исходя из физических характеристик задачи).2. Делаем один шаг метода Эйлера от точки t0 в точку t0 +h0 , т. е. находимзначение X̂(t0 + h0 ).3. Вычисляем новый шаг интегрирования по формулам (4.14)vu 2(εm + η)u .h1 = u t ẍj (t0 +h0 ) m x̃j (t0 +h0 ) j=1 4. Рассмотрим возможные случаи:если h1 = h0 : вычисляем решение в точке t0 +h0 , полагаем t0 равным t0 +h0и переходим на шаг 1.если h1 < hmin или h1 > hmax : полагаем h0 = (hmin + hmax )/2 и переходимна шаг 2.если hmin < h0 < h1 < hmax : находим h2 по формулам (4.14).Если hmin < h1 < h2 < hmax , полагаем h0 равным h1 ипереходим на шаг 3.Если h2 = h1 : вычисляем решение в точке t0 + h1 , полагаем t0 равным t0 + h1 и переходим на шаг 1.Во всех остальных случаях полагаем h0 равным h1 ,полагаем hmin равным h1 и переходим на шаг 2.если hmin < h1 < h0 < hmax , находим h2 по формулам (4.14).Если hmin < h2 < h1 < hmax , полагаем h0 равным h1 ипереходим на шаг 3.Если h2 = h1 , вычисляем решение в точке t0 + h1 , полагаем t0 равным t0 + h1 и переходим на шаг 1.Во всех остальных случаях полагаем h0 равным h1 ,полагаем hmax равным h1 и переходим на шаг 2.2104.4.
Численные примерыПокажем, как работает приведенный метод. Для этого рассмотрим несколько систем обыкновенных дифференциальных уравнений.Пример 12. Для тестирования метода было проведено его сравнение с другими методами решения обыкновенных дифференциальных уравнений для одногопростого жесткого дифференциального уравнения, которое было рассмотренов статье [117]y 0 (x) = y 2 − y 3 , y(0) = 10−4 , 0 ≤ x ≤ 20000.(4.16)Рис.
4.1. Решение жесткого дифференциального уравнения предложенным методом. Криваяв точности совпадает с графиком аналитического решения.В таблице 4.1 приведено сравнение различных методов численного интегрирования задачи Коши (4.16) на промежутке [0, 20000].211Замечание 44. Использованы следующие обозначения: ML: различные программы для решения обыкновенных дифференциальных уравнений (Matlab); KT:два метода, предложенные в статье [117] (ограниченный и неограниченный);EM: метод, предложенный в диссертации; E: метод Эйлера с постояннымшагом интегрирования; R: классический метод Розенброка.Замечание 45. Программы для метода Эйлера и предлагаемого метода написаны на языке C++ и расчет произведен в арифметике с плавающей точкойс одинарной точностью.
Правая часть вычислялась на каждом шаге интегрирования.На рис. 4.1 представлен график, полученный с помощью предлагаемогометода, который в точности совпадает с графиком аналитического решения.Пример 13. Рассмотрим задачу Коши для жесткой системы ОДУ3xdx= 106 x + y −,dt3dy= −x.dt(4.17)В таблице 4.2 приведено сравнение различных методов численного интегрирования задачи Коши (4.17) на промежутке [0, 10].МетодВремя счета, сМетодВремя счета, сEM0.1089ML ode23s0.12 ± 0.01ML ode15s0.12 ± 0.04ML ode451.04 ± 0.04ML ode23t0.13 ± 0.02KT, ограниченный2.10 ± 0.03ML ode23tb0.15 ± 0.02KT, неограниченный4.08 ± 0.01E, шаг 0.022656R, шаг 30.219Таблица 4.1.
Сравнение различных методов решения ОДУ для жесткого уравнения (4.16)212Замечание 46. Через EM обозначен предложенный в диссертации метод; E— метод Эйлера с постоянным шагом интегрирования; R — классический метод Розенброка; RK — метод Рунге — Кутта. Промежуток [0, 2] — меньшийпромежуток интегрирования для той же задачи.Замечание 47.
Программы для вычисления методом Эйлера и методом, рассмотренным в диссертации, написаны на языке C++, вычисления производились с одинарной точностью в арифметике с плавающей точкой. Правыечасти системы вычислялись на каждом шаге интегрирования.На рис. 4.2 представлен график решения, полученный методом, предложенным в диссертации.На рис.
4.3 представлены решения системы уравнений (4.17) на промежутке [0, 2]. Одно из них получено с помощью приведенного в диссертации метода,другое — с помощью метода Эйлера с постоянным шагом интегрирования.Пример 14. Рассмотрим систему обыкновенных дифференциальных уравнений, представляющую собой модель типа Ходжкина — Хаксли [64].МетодВремя счета,МетодсВремя счета,сEM3.969R, шаг 10−517.141RK, шаг 10−528.813R, шаг 2 · 10−5расходитсяEM, [0, 2]0.781E, [0, 2],расходитсяшаг 10−6E, [0, 2],шаг 10−8расходитсяE, [0, 2],1.312шаг 10−7Таблица 4.2.
Сравнение различных методов решения ОДУ для жесткой системы (4.17)213Рис. 4.2. Решение жесткой системы ОДУ (4.17) предложенным методом на промежутке [0, 10].Полученная кривая в точности совпадает с графиками наиболее точных решений, полученных другими методами.Имеем следующие уравнения модели:ṁ = (m∞ (V ) − m)/τm ,(4.18)ṅ = (n∞ (V ) − n)/τn ,(4.19)ḣ = (h∞ (V ) − h)/τh ,(4.20)IK = ḡK n4 (V − EK ),INa = ḡNa m3∞ (V )(1 − n)(V − ENa ),(4.21)(4.22)214Рис. 4.3.
Решения системы ОДУ (4.17) на промежутке [0, 2] (сплошная линия: метод, предложенный в диссертации; штриховая линия: метод Эйлера с постоянным шагом интегрирования 10−7 ).INaPh = ḡNaPh p∞ (V )h(V − ENa ),(4.23)IL = ḡL (V − EL ),(4.24)Im = INa + IK + INaPh + IL − Istim ,(4.25)V̇ = −Im /Cm .(4.26)Здесь V — потенциал мембраны, Istim — значение приложенной силы тока,Im — общая сила тока, текущего через мембрану. Потенциалозависимые установившиеся функции активации (m∞ (V ), n∞ (V ), h∞ (V ), p∞ (V )) имеют видx∞ (V ) = 1/(1 + exp((V − θx )/dx )).215Рис. 4.4. Зависимость V от времени (a: метод Эйлера с постоянным шагом интегрирования0.1 мс; b: метод, описанный в диссертации).Синхронизирующие переменные параметры приведены в Таблице 4.3.
Дополнительные параметры таковы: ḡNa = 28 нс, ḡK = 11.2 нс, ḡNaPh = 2.8 нс,ḡL = 2.8 нс, ENa = 50 мВ, EK = −85 мВ, EL = −65 мВ, Cm = 28 пФ.Переменная x θx (mV) dx (mV) τx (мс)m−34−50.1n−29−410.0h−486p−40−610000NAТаблица 4.3. Параметры для синхронизирующих переменных дифференциальных уравнений216Рис. 4.5. Зависимость Istim от времени.Начальные данные были выбраны следующими:m = 0.02, n = 0.0, h = 0.544, V = −55.0 мВ, Istim = 18 пА.Сначала оценим относительную погрешность η, т. е. норму вектора, компоненты которого равны относительным погрешностям округления, возникающимпри вычислении правых частей системы уравнений (4.8).Найдем относительную погрешность x∞ (V ). Относительная погрешностьexp(α) равна ε (см. [135]).
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
