Если у вас возникли сложности с курсовой, контрольной, дипломной, рефератом, отчетом по практике, научно-исследовательской и любой другой работой - мы готовы помочь.
Предоплата всего
Подписываем
8
Численные методы решения обыкновенных
дифференциальных уравнений.
Настоящая работа ставит целью сравнительный анализ вычислительных схем решения задачи Коши для обыкновенных дифференциальных уравнений (ОДУ). Такая задача формули руется в следующем виде
(1)
Далее везде будем предполагать, что условия существования решения задачи (1) выполнены. Также, там, где это не принципиально, будем рассматривать методы решения скалярной задачи (1).
1. Вычислительные схемы
В принципе, наиболее простым способом построения решения в точке t, если оно известно в точке t, является способ, основанный на разложении функции решения в ряд Тейлора
(2)
где
Если этот ряд оборвать и заменить приближенным значением , то получим приближенную формулу:
При р=1 она представляет собой вычислительную схему явного метода Эйлера
(3)
Применение формулы (2) ограничено лишь теми задачами, где легко вычисляется производные высших порядков функции f(y,t) правой части уравнения (1). Заметим, что обычно это не так.
В начале нашего века Рунге, Хойн и Кутта предложили подход, осно ванный на построении формулы для вида
(4)
в которой функция Ф близка к F , но не содержит производных от функции правой части уравнения. Было получено семейство явных и неявных методов, требующих s-кратного вычисления функции правой части на каждом шаге интегрирования (s-этапные методы).
Формулы этих методов идеально приспособлены для практических расчетов: они позволяют легко менять шаг интегрирования, являются одношаговыми, достаточно экономичны, по крайней мере, до формул четвертого порядка включительно. Возможно, наиболее известной является формула четырех этапного метода четвертого порядка:
(5)
Одна из основных проблем, связанных с применением методов (3), (5), (а в действительности, всех явных методов) состоит в выборе величины шага интегрирования h, обеспечивающей устойчивость вычисли тельной схемы.
За счет некоторого вычислительного усложнения этих формул был получен класс неявных методов, у которых отмеченная проблема в значительной степени снята. Неявные вычислительные схемы представляют собой алгебраические уравнения, в общем случае нелинейные, относительно значений . Например:
(6)
- метод трапеции
(7)
Получив приближение к решению в точках t, t,.., t можно использовать их для нахождения решения в точке t. Эта идея приводит к группе многошаговых методов, например, к методам Адамса. В форме Лагранжа вычислительные схемы неявных методов этого класса (методы Адамса-Мултона) имеют вид:
(8)
Заметим, что неявный метод Эйлера и метод трапеции являются частными случаями последней вычислительной схемы. Очевидно, использование многошаговых формул ставит задачу вычисле ния р штук начальных значений y , y... y , точность задания которых должна быть не хуже точности соответствующей формулы.
Отмеченная выше трудность выбора шага интегрирования h, обеспе чивающего численную устойчивость метода, делает неявные схемы предпоч тительнее для практического использования. В этой связи отметим некоторые проблемы их реализации.
Из (8) получаем:
(9)
Можно показать [1], что если выполняется неравенство
,
где L - константа Липшица для функции f(y), то существует единственное решение алгебраического уравнения (9) , которое можно получить методом простой итерации:
(10)
Очевидна желательность выбора "хорошего" начального приближения . Его надлежащий выбор обеспечивается посредством явной формулы того же порядка точности. В этом случае явная схема выполняет роль прогнозирующей, а неявная формула (10) реализует коррекцию решения и весь комбинированный процесс становится методом прогноза-коррекции - П(ВК)k, где k – номер итерации в процессе (10).
Реализация неявных вычислительных методов по схеме (10) не являет ся единственно возможной. Для решения получающихся нелинейных уравне ний и их систем (как в методах Рунге-Кутты) широко применяется метод Ньютона. Однако, проблема вычисления начального приближения и общая идея схемы «прогноз-коррекция» остается без изменений. В этой связи обра тим внимание на то, что переход к неявным схемам в значительной степени снял проблему выбо ра величины шага интегрирования h как фактора, определяющего устойчивость метода, но привел к проблемам выбора начального приближения и величины шага h, обеспечивающих сходимость итерационного процесса решения нели нейных алгебраических уравнений. Ясно, что для линейных задач этих проблем нет и применение неявных вычислительных схем для них существен но проще и надежней.
Наличие жестких задач (жесткость - свойство задачи, а не метода) сделало неявные вычислительные схемы особенно привлекательными и привело к разработке группы специальных жестко-устойчивых методов (мето ды Гира). Напомним, что явление жест кости возникает в системах ОДУ, являющихся моделями физических систем, в которых протекают процессы с существенно (на несколько порядков) от личающимися постоянными времени.
Для того, чтобы определить является ли данная задача Коши жесткой, необходимы сведения о ее поведении в окрестности частного решения y(t). В окрестности такого решения уравнение (1) можно хорошо аппрок симировать линеаризованным уравнением
(11)
где - матрица Якоби векторной функции , вычисленная в точке .
Если J(t) на некотором интервале t изменяется мало и - собственные числа матрицы Якоби, то задача (11) будет жесткой в некотором интервале , если для t T выполняется:
Систему уравнений можно считать жесткой, если больше 10, однако во многих задачах радиоэлектроники, физики, химической кинетики, управления этот коэффициент достигает значений порядка 106 и более.
Введение понятия жесткой устойчивости метода (Гир, 1969, 1971) позволило сконструировать группу методов, в которых величина шага интегрирования выбирается так, чтобы быстро затуха ющие и не оказывающие существенного влияния компоненты решения аппроксимировались устойчиво, тогда как для компонент с большими постоянными времени гарантировалась точность аппроксимации.
Наиболее распространенным классом линейных многошаговых методов для жестких задач являются "формулы дифференцирования назад" (более общее название - методы Гира). Вычислительная формула их основывается на интерполяционной формуле Лагранжа для правой части уравнения (1)
, (12)
Выражение (12) представляет собой нелинейное уравнение относитель но ym+1. Для вычисления начального приближения по значениям строится интерполяционный полином k-го порядка для функции y(t) и посредством экстраполирования вычисляется значение . Заметим, что по формулам (12) можно вести интегрирование с переменным шагом. Вместе с тем, при фиксированном шаге коэффициенты не изменяются, что существенно упрощает метод.
2. Погрешность численного решения и методы ее оценки
Главным источником погрешности численного решения в точке
(13)
является погрешность аппроксимации метода, которая возникает на каждом шаге интегрирования. Для этой погрешности, в предположении гладкости решения, справедлива оценка
, (14)
то для полной погрешности (13) справедлива оценка
(15)
где: h = max {h} , 1 i n и L – константа Липшица.
Очевидно, что даже упрощенные оценки погрешности по приведенным формулам не представляют практического интереса, поскольку требуют вычисления или оценивания верхних границ частных производных высших порядков от f(t,y) (постоянные C и L). Для практики однако, оценки погрешности необходимы, с одной стороны, чтобы реализовать выбор величины hn , достаточно малой для получения требуемой точности, и одновременно - достаточно большой для минимизации вычислительных затрат.
Самый старый способ оценки погрешности для одношаговых методов, который использовал еще Рунге, состоит в двукратном вычислении значения yn с шагом h и h/2. При этом, погрешность решения, полученного с шагом h, оценивается выражением
, (16)
а решения, полученного с шагом h/2, -
. (17)
Еще одна идея оценки величины ошибки численного решения состоит в выполнении шага от точки tn-1 к tn методами порядка р и р+1 для получения оценки ошибки интегрирования методом порядка р .
(18)
Очевидна вычислительная «дороговизна» такого способа оценки ошибки. Эта идея наиболее элегантное свое воплощение нашла в предложенных Инглендом (1967) и Фельбергом (1969) вложенных методах. Они представляют собой формулы Рунге-Кутты, которые одновременно дают возможность получить решения порядка p и р+1.
Первые методы такого типа были предложены Мерсоном (1957). Однако, в его методе решение y(p+1) имеет пятый порядок точности только для линейных уравнений с постоянными коэффициентами. Поэтому этот метод в общем случае переоценивает погрешность при малых h, но все же работает вполне удовлетворительно.
Формулы Фельберга не имеют подобных ограничений. Он получил целое семейство таких методов различного порядка, но наиболее популярным стал шестиэтапный метод 4-го порядка. Следует иметь в виду, что Фельберг минимизировал коэффициент погрешности аппроксимации для решения низшего порядка (4-го). Поэтому точность результата высшего (5-го) порядка трудно поддается оценке и этот результат используется только для регулирования величины шага.
Использование процедуры прогноза-коррекции для методов Гира дает простой и эффективный способ оценки погрешности аппроксимации. Брайтон (1972) показал, что ее главную часть можно выразить через разность предсказанного и скорректированного значений решения. Так, для метода k-го порядка справедлива оценка
(19)
3. Общие проблемы реализации численных методов
Одной из центральных проблем реализации алгоритмов численного интегрирования задачи Коши является выбор стратегии регулирования шага интегрирования и оценки допустимой его величины при интегрировании с постоянным шагом.
При интегрировании с постоянным шагом выбор величины шага определяется, главным образом, условиями устойчивости метода. Для линейных задач исчерпывающую информацию на этот счет несут собственные числа матрицы Якоби правой части системы. Для нелинейных задач максимальный шаг можно оценить величиной , где L – константа Липшица функции и С – постоянная, существенно зависящая от метода, но редко превосходящая 10.
Методы, обладающие более высокими показателями устойчивости являются неявными и их реализация приводит к необходимости решения систем нелинейных алгебраических уравнений. Ограничения, накладываемые свойствами устойчивости метода, при этом, могут быть существенно ужесточены условиями сходимости итерационного процесса, применяемого для решения алгебраических систем уравнений, порождаемых численным методом интегрирования. Метод простой итерация вида (10) часто практически неприемлем (для жестких задач особенно), так как для его сходимости требуется выполнение условия, накладывающего на величину шага столь же сильное ограничение, как и то, которого старались избежать, переходя от явных схем к неявным.
В таких ситуациях для решения нелинейных уравнений оказывается более эффективным метод Ньютона, который позволяет нередко добиться сходимости при существенно меньших ограничениях на величину шага. Точность прогноза для , используемого для начала итерационного ньютоновского процесса, может оказывать заметное влияние на скорость сходимости. Поэтому следует иметь в виду, что дисперсия величины прогноза быстро растет c ростом числа значений i = 1,2,…,k, участвующих в предсказывающей формуле, что может привести к существенному снижению эффективности многошаговых формул при k > 4.
Если предпринять все необходимые меры для снятия ограничения на величину шага по условиям устойчивости метода и сходимости итерационного процесса, то выбор шага может быть обусловлен лишь требуемой точностью решения. При решении жестких систем возможны две стратегии выбора шага.
Для эффективного использования последней стратегии необходим относительно простой механизм изменения шага интегрирования в соответствии с оценками погрешности аппроксимации, приведенными выше. На основании формул (16 - 19) вычисляется оценка ошибки решения
где - масштабирующий множитель для i-й компоненты вектора решения (при вычислении абсолютной погрешности - =1, а при вычислении относительной погрешности - ). Для повышения надежности программы разумно использовать масштабирование вида
Величина err сравнивается с заданной допустимой погрешностью Tol и на основании оценки (14) вычисляется оптимальная величина следующего шага
Соблюдая определенную осторожность, этот шаг следует несколько уменьшить умножением на 0.8 или 0.9. Не следует допускать и слишком быстрых изменений величины шага – разумно ограничить его вариации величиной .
Если , то решение, полученное с шагом , следует считать удовлетворительным и процесс интегрирования можно продолжить с шагом . В противном случае решение с шагом отбрасывается, и интегрирование из точки начинается с шагом .
Значительные трудности возникают при определении величины Tol – допустимой ошибки на одном шаге. Брайтон с соавторами предложил определять ее из условия
где - предельная допустимая величина ошибки на полном интервале интегрирования T. Заметим, что такая стратегия определения этой величины для очень многих, особенно нежестких, задач приводит к чрезмерно сильным ограничениям величины шага.
4. Описание программного комплекса COSHI.
Программа COSHI позволяет провести изучение характеристик методов численного интегрирования систем ОДУ в форме Коши. Она включает в себя следующие методы:
Неявные методы Эйлера, трапеции и Гира реализованы в общем виде и в вариантах, предназначенных для интегрирования систем линейных уравнений с постоянным шагом. Такие ограничения позволяют построить достаточно эффективную процедуру, суть которой продемонстрирована ниже на примере метода трапеции.
Линейная задача Коши представима в форме:
Стандартная вычислительная схема метода трапеции в этом случае легко преобразуется в выражение
При постоянном шаге h эта схема требует лишь однократного формирования матриц (Е-0.5hА) и (Е+0.5hА) и лишь однократного обращения матрицы (Е-0.5hА). Таким образом, вычисление вектора решения во всех точках интервала интегрирования сводится к достаточно экономичным операциям умножения и суммирования векторов и матриц.
Реализация стандартной схемы метода трапеции представлена двумя процедурами: в одной на этапе коррекции использована простая итерационная схема П(ВК)k, в другой – схема ньютоновских итераций до сходимости. В последней процедуре предусмотрена возможность регулирования шага по методу повторного интегрирования с половинным шагом (выражение (16).
Методы Рунге-Кутты запрограммированы в варианте без регулирования шага.
Алгоритмы Гира представлены в модификациях для линейных систем с постоянным шагом интегрирования и в общем виде, позволяющем регулировать шаг по схеме (19). При этом, регулирование шага осуществляется в соответствии с заданной точностью решения (). При этом допустимая погрешность на шаге интегрирования устанавливается в соответствии с рекомендацией Брайтона.
В программе установлена верхняя граница для величины шага интегрирования . Предприняты также меры, предупреждающие быстрое изменение величины шага. Все это довольно часто приводит к существенно меньшей результирующей погрешности интегрирования, чем задается пользователем, но повышает надежность работы программного комплекса.
Программа позволяет записать решение для любых трех компонент вектора и представить их графически. Существует возможность изображения фазового портрета процесса относительно двух, из указанных выше трех, компонент решения. При решении линейных задач вычисляются собственные числа матрицы задачи, что позволяет правильно выбрать шаг интегрирования и спрогнозировать характер решения.
Исследование свойств методов может быть проведено на тестовых задачах, включенных в программный комплекс, и на задачах пользователя. Последние должны быть включены в пакет посредством описания их в процедурах bUse, VectFunUs и MatUse, находящихся в файле UsesFl.pas.
5. Тестовые задачи
3. Жесткая задача.
4. Задача с регулируемой жесткостью
5.
6.
7.
11. Нелинейная задача с регулируемой жесткостью.
12.
13.
14.
15.
16.
Задачи с порядковым номером 21 и 31 записаны в процедурах MatrUses, bUse и VectFunUs (файла UsesFl) в качестве примера оформления пользовательских линейных и нелинейных задач.
21.
31.
Под номерами 22 и 32 могут записываться линейные и нелинейные (соответственно) пользовательские задачи. Включение этих задач производится в среде TURBO-PASCAL и предполагает перекомпиляцию программного пакета. При этом пользователь работает только с файлом UsesFl.Pas, указав в директиве главного меню «Compile» в опции Primeri File имя «Coshi.pas». Компиляция производится по директиве «Make».
6. Задание
1. При интегрировании нескольких (5-6) систем линейных уравнений на основе информации о значениях собственных чисел матрицы системы оценить максимальную величину шага устойчивого интегрирования и проверить эту оценку экспериментально. Зафиксировать величину шага, при которой метод теряет устойчивость.
2. При интегрировании нескольких (3-4) систем нелинейных уравнений обратить внимание на эффективность различных методов реализации схемы прогноз-коррекция.
3. Исследовать поведение полной ошибки численного решения при интегрировании с постоянным шагом методами различного порядка. (Методы Эйлера, РК2 – РК4), т.е. получить зависимость максимальной погрешности решения задачи от порядка метода при условии постоянства шага интегрирования, и методами одного порядка при вариациях величины шага. Дать качественное описание поведения функции полной погрешности решения.
Содержание отчета.
Отчет должен включать: