Математический форум Math Help Planet
Обсуждение и решение задач по математике, физике, химии, экономике Теоретический раздел |
Часовой пояс: UTC + 3 часа [ Летнее время ] |
новый онлайн-сервис число, сумма и дата прописью |
|
Часовой пояс: UTC + 3 часа [ Летнее время ] |
Разностные схемы для решения задачи Коши | |
---|---|
Онлайн-сервисы
Нахождение НОД и НОК
Разложение числа на простые множители
Сравнения по модулю
Операции над множествами
Операции над векторами
Разложение вектора по базису. Доказательство, что векторы образуют базис
Чертёж треугольника по координатам вершин
Решение треугольника
Решение Пирамиды
Построение Пирамиды по координатам вершин
Чертёж многоугольника по координатам вершин
Решение систем методом Крамера и Матричным
Онлайн построение графика кривой 2-го порядка
Определение вида кривой или поверхности 2-го порядка по инвариантам
МНК и регрессионный анализ Онлайн + графики
Онлайн число, сумма и дата прописью
Алгоритмы JavaScript
Алгоритмы поиска
Алгоритмы сортировки
Уникальные элементы массива
Объединение, пересечение и разность массивов
НОД и НОК
Операции над матрицами
Дата прописью
Введение в анализ
Функции: понятие, определение, графики
Непрерывность функции
Исследование функции и построение графика
Теория множеств
Множества: понятие, определение, примеры
Точечные множества
Замкнутые и открытые множества
Мера множества
Группы, кольца, поля в математике
Поле комплексных чисел
Кольцо многочленов
Основная теорема алгебры и ее следствия
Математическая логика
Алгебра высказываний
Аксиоматика и логические рассуждения
Методы доказательств теорем
Алгебра высказываний и операции над ними
Формулы алгебры высказываний
Тавтологии алгебры высказываний
Логическая равносильность формул
Нормальные формы для формул высказываний
Логическое следование формул
Приложение алгебры высказываний для теорем
Дедуктивные и индуктивные умозаключения
Решение логических задач
Принцип полной дизъюнкции
Булевы функции
Множества, отношения и функции в логике
Булевы функции от одного и двух аргументов
Булевы функции от n аргументов
Системы булевых функций
Применение булевых функций к релейно-контактным схемам
Релейно-контактные схемы в ЭВМ
Практическое применение булевых функций
Теория формального
Формализованное исчисление высказываний
Полнота и другие свойства формализованного исчисления высказываний
Независимость системы аксиом формализованного исчисления высказываний
Логика предикатов
Логика предикатов
Логические операции над предикатами
Кванторные операции над предикатами
Формулы логики предикатов
Тавтологии логики предикатов
Преобразования формул и следование их предикатов
Проблемы разрешения для общезначимости и выполнимости формул
Применение логики предикатов в математике
Строение математических теорем
Аристотелева силлогистика и методы рассуждений
Принцип полной дизъюнкции в предикатной форме
Метод полной математической индукции
Необходимые и достаточные условия
Логика предикатов и алгебра множеств
Формализованное исчисление предикатов
Неформальные и формаль-ные аксиоматические теории
Неформальные аксиоматические теории
Свойства аксиоматических теорий
Формальные аксиоматические теории
Формализация теории аристотелевых силлогизмов
Свойства формализованного исчисления предикатов
Формальные теории первого порядка
Формализация математической теории
Теория алгоритмов
Интуитивное представление об алгоритмах
Машины Тьюринга и тезис
Рекурсивные функции
Нормальные алгоритмы Маркова
Разрешимость и перечислимость множеств
Неразрешимые алгоритмические проблемы
Теорема Гёделя о неполноте формальной арифметики
Математическая логика и компьютеры
Дискретная математика
Множества и отношения
Теория множеств: понятия и определения
Операции над множествами
Кортеж и декартово произведение множеств
Соответствия и бинарные отношения на множествах
Операции над соответствиями на множествах
Семейства множеств
Специальные свойства бинарных отношений
Отношения эквивалентности на множестве
Упорядоченные множества
Теорема о неподвижной точке
Мощность множества
Парадокс Рассела
Метод характеристических функций
Группы и кольца
Алгебраические структуры и операции
Группоиды, полугруппы, группы
Кольца, тела, поля
Области целостности в теории колец
Модули и линейные пространства
Подгруппы и подкольца
Теорема Лагранжа о порядке конечной группы
Гомоморфизмы групп и нормальные делители
Гомоморфизмы и изоморфизмы колец
Алгебра кватернионов
Полукольца и булевы алгебры
Полукольца: определение, аксиомы, примеры
Замкнутые полукольца
Полукольца и системы линейных уравнений
Булевы алгебры и полукольца
Решетки и полурешетки
Алгебраические системы
Алгебраические системы: модели и алгебры
Подсистемы алгебраических систем
Конгруэнции и фактор-системы
Гомоморфизмы алгебраических систем
Прямые произведения алгебраических систем
Конечные булевы алгебры
Многосортные алгебры
Теория графов
Теория графов: основные понятия и определения
Способы представления графов
Неориентированные и ориентированные деревья
Остовное дерево и алгоритм Краскала
Методы систематического обхода вершин графа
Алгоритмы поиска в глубину и ширину в графах
Задача о путях во взвешенных ориентированных графах
Изоморфизм, гомоморфизм и автоморфизм графов
Топологическая сортировка вершин графа
Элементы цикломатики в теории графов
Булева алгебра и функции
Булевы функции и булев куб
Таблицы булевых функций и булев оператор
Равенство булевых функций. Фиктивные переменные
Формулы и суперпозиции булевых функций
Дизъюнктивные и конъюнктивные нормальные формы
Построение минимальных ДНФ
Теорема Поста и классы
Критерий Поста
Схемы из функциональных элементов
Конечные автоматы и регулярные языки
Конечные автоматы и регулярные языки
Алфавит, слово, язык в программировании
Порождающие грамматики (грамматики Хомского)
Классификация грамматик и языков
Регулярные языки и регулярные выражения
Конечные автоматы
Допустимость языка конечным автоматом
Теорема Клини
Детерминизация конечных автоматов
Минимизация конечных автоматов
Лемма о разрастании для регулярных языков
Обоснование алгоритма детерминизации автоматов
Конечные автоматы с выходом
Морфизмы и конечные подстановки
Машины Тьюринга
Контекстно-свободные языки
Контекстно-свободные языки и грамматики
Приведенная форма КС-грамматики
Лемма о разрастании для КС-языков
Магазинные автоматы (автомат с магазинной памятью)
Алгоритм построения МП-автомата по КС-грамматике
Алгоритм построения КС-грамматики по МП-автомату
Алгебраические свойства КС-языков
Основное свойство суперпозиции КС-языков
Пересечение контекстно-свободных языков
Методы синтаксического анализа КС-языков
Восходящий синтаксический анализ и LR(k)-грамматики
Семантика формальных языков
Принцип индукции по неподвижной точке
Графовое представление МП-автоматов
Интегральное исчисление
Неопределённый и определённый
Неопределенный и определенный интегралы
Свойства интегралов
Интегрирование по частям
Интегрирование методом замены переменной
Интегрирование различных рациональных функций
Интегрирование различных иррациональных функций
Интегрирование различных тригонометрических функций
Определенный интеграл и его основные свойства
Необходимое и достаточное условие интегрируемости
Теоремы существования первообразной
Свойства определенных интегралов
Несобственные интегралы
Интегральное определение логарифмической функции
Приложения интегралов
Вычисление площадей плоских фигур
Площади фигур в различных координатах
Вычисление объемов тел с помощью интегралов
Объём тела вращения
Вычисление длин дуг кривых
Формулы длины дуги регулярной кривой
Кривизна плоской кривой
Площадь поверхности вращения тела
Интегралы в физике
Статические моменты и координаты центра тяжести
Теоремы Гульдина–Паппа
Вычисление моментов инерции
Другие приложения интегралов в физике
Основные интегралы
Вариационное исчисление
Примеры вариационных задач
Дифференциальное уравнение Эйлера
Функционалы, зависящие от нескольких функций
Задача о минимуме кратного интеграла
Финансовый анализ
Анализ эффективности
Критерии и показатели эффективности предприятия
Методы анализа эффективности деятельности
Факторный анализ прибыли от операционной деятельности
Анализ безубыточности предприятия
Операционный рычаг и эффект финансового рычага
Анализ и оценка состава, структуры и динамики доходов и расходов
Анализ рентабельности и резервов устойчивого роста капитала
Анализ распределения прибыли предприятия
Анализ и оценка чувствительности показателей эффективности
Анализ устойчивости
Финансовая устойчивость и долгосрочная платежеспособность
Характеристика типов финансовой устойчивости
Рыночная активность
Финансовый анализ рыночной активности
Методика анализа рыночной активности
Анализ и оценка дивидендного дохода на одну акцию
Инвестиционная деятельность
Инвестиции: экономическая сущность и классификация
Государственное регулирование инвестиционной деятельности
Источники финансовых ресурсов на капитальные вложения
Инвестиции в основные фонды
Оценка состояния основных фондов
Амортизация основных фондов
Капитальное строительство в инвестиционном процессе
Планирование инвестиций в форме капитальных вложений
Экономическая эффективность инвестиций
Финансирование капитальных вложений
Кредитование капитальных вложений
Кредитоспособность
Финансирование и кредитование затрат
Финансирование и кредитование инвестиционной деятельности потребительской кооперации
Финансирование и кредитование капитальных вложений потребительской кооперации
Инвестиционное строительное проектирование
Анализ инвестиций
Инвестиции и инвестиционная деятельность предприятия
Задачи финансового анализа инвестиций предприятия
Учет фактора времени в инвестиционной деятельности
Аннуитет и финансовая рента в инвестициях
Учет фактора инфляции при инвестировании
Оценка фактора риска инвестиционного проекта
Методы оценки эффективности инвестиций
Показатели эффективности инвестиционного проекта
Стоимость компании
Концепция построения международных стандартов финансовой отчетности (МСФО)
Экономическое содержание международных стандартов финансовой отчётности
Цели и принципы оценки стоимости акций и активов компании
Оценка акций и активов предприятия по справедливой стоимости
Методы оценки справедливой стоимости акций предприятия
Затратный подход к оценки стоимости компаний и акций
Сравнительный подход к оценки стоимости предприятий и акций
Доходный подход к оценке стоимости компании и акций
Выбор ставки дисконтирования при инвестировании в акции
Метод капитализации прибыли
Сравнение подходов к оценке стоимости компаний и пакетов акций
Форвардные контракты
Форвардный контракт и цена
Форвардная цена акции на бирже
Цена форвардного контракта инвестора
Форвардная цена акции с учетом величины дивиденда
Форвардная цена акции с учетом ставки дивиденда
Форвардная цена валюты на рынке форекс
Форвардный валютный курс и инфляция на рынке
Форвардная цена товара и спотовый рынок
Форвардная цена при различии ставок по кредитам и депозитам
Синтетический форвардный контракт на акции и валюту
Теория вероятностей
Основные понятия теории вероятностей
Зависимые и независимые случайные события
Повторные независимые испытания
Формула Бернулли
Одномерные случайные величины
Многомерные случайные величины
Функции случайных величин
Законы распределения целочисленных случайных величин
Законы распределения непрерывных случайных величин
Предельные теоремы теории вероятностей
Закон больших чисел и предельные теоремы
Вероятностные закономерности
Математическая статистика
Элементы математической статистики
Выборочный метод
Оценки параметров генеральной совокупности
Статистические гипотезы
Критерии согласия
Теоретические и эмпирические частоты
Теория очередей (СМО)
Определение системы массового обслуживания
Уравнения Колмогорова
Предельные вероятности состояний
Определение СМО с отказами
Определение СМО с ожиданием (очередью)
Аналитическая геометрия
Векторная алгебра
Метрические понятия и аксиомы геометрии
Равенство и подобие геометрических фигур
Бинарные отношения
Вектор, его направление и длина
Линейные операции над векторами
Линейная зависимость и независимость векторов
Отношение коллинеарных векторов
Проекции векторов на прямую и на плоскость
Угол между векторами
Ортогональные проекции векторов
Координата вектора на прямой и базис
Координаты вектора на плоскости и базис
Координаты вектора в пространстве и базис
Операции над векторами в координатной форме
Ортогональный и ортонормированный базисы
Cкалярное произведение векторов и его свойства
Выражение скалярного произведения через координаты векторов
Векторное произведение векторов и его свойства
Смешанное произведение векторов и его свойства
Ориентированные площади и объемы
Двойное векторное произведение и его свойства
Применение векторов в задачах на аффинные свойства фигур
Применение произведений векторов при решении геометрических задач
Применение векторной алгебры в механике
Системы координат
Прямоугольные координаты
Преобразования прямоугольных координат
Полярная система координат
Цилиндрическая система координат
Сферические координаты
Аффинные координаты
Аффинные преобразования координат
Аффинные преобразования плоскости
Примеры аффинных преобразований плоскости
Аффинные преобразования пространства
Многомерное координатное пространство
Линейные и аффинные подпространства
Скалярное произведение n-мерных векторов
Преобразования систем координат
Геометрия на плоскости
Алгебраические линии на плоскости
Общие уравнения геометрических мест точек
Алгебраические уравнения линий на плоскости
Уравнения прямой, проходящей через точку перпендикулярно вектору
Уравнения прямой, проходящей через точку коллинеарно вектору
Уравнения прямой, проходящей через две точки
Уравнения прямой с угловым коэффициентом
Взаимное расположение прямых
Примеры задач с прямыми на плоскости
Системы неравенств с двумя неизвестными
Системы линейных уравнений с двумя неизвестными
Линии 2-го порядка
Канонические уравнения линий второго порядка
Порядок приведения уравнения линии к каноническому виду
Эллипс
Гипербола
Парабола
Квадратичные неравенства с двумя неизвестными
Применение линий 1-го и 2-го порядков в задачах на экстремум функций
Инварианты линий
Классификация линий 2-го порядка по инвариантам
Приведение уравнения линии к каноническому виду по инвариантам
Геометрия в пространстве
Способы задания ГМТ в пространстве
Алгебраические уравнения поверхностей
Уравнения плоскости, проходящей через точку перпендикулярно вектору
Уравнения плоскости, компланарной двум неколлинеарным векторам
Уравнения плоскости, проходящей через три точки
Взаимное расположение плоскостей
Типовые задачи с плоскостями
Уравнения прямых в пространстве
Взаимное расположение прямых в пространстве
Типовые задачи с прямыми в пространстве
Поверхности 2-го порядка
Канонические уравнения поверхностей
Порядок приведения уравнения поверхности к каноническому виду
Поверхности второго порядка
Эллипсоиды
Гиперболоиды
Конусы
Параболоиды
Применение поверхностей 1-го и 2-го порядков в задачах на экстремум функций
Инварианты поверхностей
Линейная алгебра
Матрицы и операции
Линейные операции над матрицами
Умножение матриц
Возведение матриц в степень
Многочлены от матриц
Транспонирование и сопряжение матриц
Блочные матрицы
Произведение и сумма матриц Кронекера
Метод Гаусса приведения матрицы к ступенчатому виду
Элементарные преобразования матриц
Определители
Определители матриц и их основные свойства
Формула полного разложения определителя
Формула Лапласа полного разложения определителя
Определитель произведения матриц
Методы вычисления определителей
Ранг матрицы
Линейная зависимость и линейная независимость строк (столбцов) матрицы
Ранг матрицы и базисный минор матрицы
Методы вычисления ранга матрицы
Ранг системы столбцов (строк)
Обратная матрица
Обратные матрицы и их свойства
Ортогональные и унитарные матрицы
Способы нахождения обратной матрицы
Матричные уравнения
Односторонние обратные матрицы
Скелетное разложение матрицы
Полуобратная матрица
Псевдообратная матрица
Системы уравнений
Системы линейных алгебраических уравнений
Метод Гаусса решения систем линейных уравнений
Структура общего решения системы уравнений
Решение систем с помощью полуобратных матриц
Псевдорешения системы линейных уравнений
Функциональные матрицы
Функциональные матрицы скалярного аргумента
Производные матриц по векторному аргументу
Линейные и квадратичные формы и их преобразования
Приведение форм к каноническому виду
Закон инерции вещественных квадратичных форм
Знакоопределенность форм вещественных квадратичных
Формы и исследование функций на экстремум
Многочленные матрицы
Многочленные матрицы (лямбда-матрицы)
Операции над лямбда-матрицами
Простые преобразования многочленных матриц
Инвариантные множители многочленной матрицы
Функции от матриц
Собственные векторы и значения матрицы
Подобие числовых матриц
Характеристический многочлен матрицы
Минимальный многочлен матрицы
Теорема Гамильтона-Кэли
Жорданова форма матрицы
Приведение матрицы к жордановой форме
Многочлены от матриц
Применение многочленов от матриц
Функции от матриц
Линейные пространства
Линейные пространства: определение и примеры
Линейная зависимость и независимость n-мерных векторов
Размерность и базис линейного пространства
Преобразования координат в линейном пространстве
Изоморфизм линейных пространств
Подпространства
Подпространства линейного пространства
Пересечение и сумма подпространств
Способы описания подпространств
Нахождение дополнения и суммы подпространств
Нахождение пересечения подпространств
Линейные отображения
Линейные многообразия
Линейные отображения
Матрица линейного отображения
Ядро и образ линейного отображения
Линейные операторы
Линейные операторы (преобразования)
Инвариантные подпространства
Собственные векторы и значения оператора
Свойства собственных векторов операторов
Канонический вид линейного оператора
Методика приведения линейного преобразования к каноническому виду
Евклидовы пространства
Евклидовы пространства
Ортогональные векторы евклидова пространства
Ортогональный базис евклидова пространства
Ортонормированный базис евклидова пространства
Ортогональные дополнения в евклидовом пространстве
Задача о перпендикуляре
Матрица и определитель Грама и его свойства
Линейные преобразования евклидовых пространств
Канонический вид ортогонального оператора евклидова пространства
Сопряженные операторы евклидова пространства
Самосопряженные операторы евклидова пространства
Приведение квадратичной формы к главным осям
Унитарные пространства и их линейные преобразования
Комплексный анализ
Комплексные числа
Комплексные числа в алгебраической форме
Комплексные числа в тригонометрической и показательной формах
Множества на комплексной плоскости
Последовательности и ряды комплексных чисел
Комплексные функции
Функции комплексного переменного. Предел, непрерывность и производная
Элементарные функции комплексного переменного
Дифференцирование функций комплексного переменного
Аналитические функции и их свойства
Конформные отображения
Функциональные ряды в комплексной области
и их свойства Интегрирование функций комплексного переменного
Функциональные ряды и последовательности
Степенные ряды и их свойства
Разложение функций в степенные ряды
Нули аналитических функций
Ряд Лорана и разложение функций по целым степеням
Особые точки, Вычеты
Изолированные особые точки функций и полюсы
Вычеты и их применение
Вычисление интегралов с помощью вычетов
Вычеты и расположение нулей многочлена
Операционное исчисление
Дифференциальные уравнения
ДУ первого порядка
Основные понятия и определения ДУ
Метод изоклин для ДУ 1-го порядка
Метод последовательных приближений
ДУ с разделяющимися переменными
Однородные ДУ
Линейные ДУ 1-го порядка
Дифференциальное уравнение Бернулли
ДУ в полных дифференциалах
Интегрирующий множитель
ДУ, не разрешенные относительно производной
Дифференциальное уравнение Риккати
Составление ДУ семейств линий
Задачи на траектории
Особые решения ДУ
ДУ высших порядков
Понятия и определения ДУ высших порядков
ДУ, допускающие понижение порядка
Линейная независимость функций
Определители Вронского и Грама
Однородные и неоднородные дифференциальные уравнения
Задача Коши и Уравнение Эйлера
Линейные ДУ с переменными коэффициентами
Метод Лагранжа решения ДУ
Краевые задачи для ДУ высших порядков
Разложение решения ДУ в степенной ряд
Разложение решения ДУ в обобщенный степенной ряд
Нахождение периодических решений ДУ
Асимптотическое интегрирование ДУ
Системы ДУ
Системы ДУ: понятия и определения
Сведение системы ДУ к одному уравнению
Нахождение интегрируемых комбинаций
Интегрирование однородных линейных систем ДУ
Методы интегрирования неоднородных систем ДУ
Преобразование Лапласа и решение ДУ и систем
Теория устойчивости
Численные методы
Методы алгебры
Численные методы линейной алгебры
Численные методы решения СЛАУ
Итерационный метод Шульца обратной матрицы
Методы решения задач о собственных значениях и векторах матрицы
Методы решения нелинейных уравнений
Методы решения систем нелинейных уравнений
Методы теории приближений
Методы приближения сеточных функций
Методы функциональной интерполяции
Методы интегрально-дифференциальной интерполяции
Методы интегрального сглаживания
Методы интерполяции и сглаживания сплайнами
Методы численного дифференцирования и интегрирования
Методы численного дифференцирования
Методы численного интегрирования
Методы решения обыкновенных ДУ
Численные методы решения задачи Коши
Разностные схемы для решения задачи Коши
Составные схемы для решения задачи Коши
Экстраполяционные методы решения задачи Коши
Непрерывно-дискретные методы решения задачи Коши
Численные методы решения краевых задач
Методы решения ДУ в частных производных
Численные методы решения уравнений математической физики с двумя переменными
Принципы построения разностных схем для уравнений в частных производных
Разностные схемы решения уравнений в частных производных 1-го порядка
Разностные схемы решения уравнений в частных производных 2-го порядка
Численные методы решения уравнений в частных производных
Численные методы решения уравнений математической физики с тремя переменными
|
Разностные схемы для решения задачи КошиРассмотрим основные принципы конструирования разностных схем. 1. Производная [math]\left.{\left(\frac{dy}{dx}\right)\!}\right|_{x=x_{i}}[/math] аппроксимируется на заранее выбранном шаблоне [math](x_{i-k+1}, x_{i-k+2},\ldots,x_{i+1})[/math], а правая часть [math]f(x,y)[/math] обыкновенного дифференциального уравнения заменяется сеточным представлением [math]f(x_{i},\widehat{y}_{i})[/math] в точке [math](x_{i},\widehat{y}_{i})[/math]. Условно этот принцип называется аппроксимационным. С его помощью строятся одношаговые и многошаговые схемы в зависимости от шаблона, по которому аппроксимируется производная. 2. Дифференциальное уравнение заменяется разностным аналогом на одношаговом (двухточечном) шаблоне [math](x_{i},x_{i+1})[/math] (возможно использование дополнительных точек внутри шаблона) путем согласования с разложением функции [math]y=y(x)[/math] при [math]x=x_{i+1}[/math] no формуле Тейлора относительно точки [math]x=x_{i}[/math]. По этому принципу строятся одношаговые схемы, в частности семейство схем Рунге-Кутты. 3. Осуществляется интегрирование обеих частей дифференциального уравнения вдоль решения на выбранном шаблоне. При этом правая часть уравнения предварительно заменяется интерполяционным многочленом, степень которого должна соответствовать порядку точности схемы. Данный принцип называется интегрально-интерполяционным и используется для построения одношаговых и многошаговых явных и неявных схем, в том числе схем Адамса. 4. Разностные схемы строятся на основе подобия соотношений теории приближения конструкциям схем. С этой целью используются следствия из интегрально-дифференциальных сплайнов — параметрические соотношения, формулы аппроксимации производных, квадратурные формулы, записанные в общем случае на нерегулярном шаблоне. Данный принцип, который можно назвать принципом аналогий, далее будет применен для конструирования двух- и трехшаговых явных и неявных схем. Принцип аппроксимацийРассмотрим проблему нахождения численного решения задачи Коши (6.9): [math]\frac{dy}{dx}= f(x,y),~ y(x_0)=y_0[/math]. Для аппроксимации производной [math]\left.{\left(\frac{dy}{dx}\right)\!}\right|_{x=x_{i}}[/math] используем формулу (5.4), записанную на двухточечном шаблоне [math]H_{2,i}= (x_{i},x_{i+1})\colon[/math] [math]\left.{\left(\frac{dy}{dx}\right)\!}\right|_{x=x_{i}}= \frac{y_{i+1}-y_{i}}{h_{i+1}}+ O(h_{i+1})\quad \left(\frac{h_{i+1}}{2}M_{2,i}\right)\!,[/math] (6.16) и формулу (5.10), записанную на трехточечном нерегулярном шаблоне [math]H_{3,i}= (x_{i-1}, x_{i}, x_{i+1})\colon[/math] [math]\begin{gathered}\left.{\left(\frac{dy}{dx}\right)\!}\right|_{x=x_{i}}= \frac{\Pi_{i}^{i+1}}{H_{i}^{i+1}}\cdot\! \left[\frac{\Delta y_{i-1}}{h_{i}^2}+ \frac{\Delta y_{i}}{h_{i+1}^2}\right]=\\ =\frac{1}{H_{i}^{i+1}}\! \left[-\delta_{i+1}y_{i-1}+ \frac{\delta_{i+1}^2-1}{\delta_{i+1}}y_{i}+ \frac{y_{i+1}}{\delta_{i+1}}\right]+ O(h_{i}h_{i+1})\quad \left(\frac{h_{i}^2}{6} \delta_{i+1}M_{3,i}\right)\!, \end{gathered}[/math] (6.17) где [math]H_{i}^{i+1}= h_{i}+h_{i+1},~ \Pi_{i}^{i+1}= h_{i}h_{i+1},~ h_{i+1}=x_{i+1}-x_{i},~~ \delta_{i+1}= \frac{h_{i+1}}{h_{i}}[/math] — параметр нерегулярности. При [math]\delta_{i+1}=1[/math] шаблон становится регулярным, а среднее слагаемое с входящим в него [math]y_{i}[/math] — равным нулю. При этом формула (6.17) сводится к традиционной: [math]\left.{\left(\frac{dy}{dx}\right)\!}\right|_{x=x_{i}}= \frac{y_{i+1}-y_{i-1}}{2h}+ O(h^2)\quad \left(\frac{h^2}{6}M_{3,i}\right)\!.[/math] (6.18) Аппроксимация формульных функций, входящих в правую часть (6.9), осуществляется путем перехода к их сеточному представлению. Методика построения разностных схем с помощью аппроксимаций производной1. Вводится в общем случае неравномерная сетка [math]\Omega_n= (x_0, x_1,\ldots, x_{i}, x_{i+1},\ldots, x_{n})[/math]. Величина шага [math]h_{i+1}=x_{i+1}-x_{i}[/math] выражается через узловые точки. 2. Для получения рассматриваемых здесь двух схем (одношаговой и двух-шаговой) выбирается соответственно двухточечный (одношаговый) и трехточечный (двухшаговый) шаблоны [math]H_{2,i}= (x_{i},x_{i+1})[/math] и [math]H_{3,i}= (x_{i-1},x_{i},x_{i+1})[/math], на которых в точке [math]x_{i}[/math] производится аппроксимация производной [math]\left.{\left(\frac{dy}{dx} \right)\!}\right|_{x=x_{i}}[/math] по формулам. (6.16), (6.17) соответственно в левой и центральной точках [math]x_{i}[/math]. Далее заменяется правая часть уравнения (6.9) ее сеточным представлением, т.е. [math]f(x,y)\to f(x_{i},y_{i})[/math], а вместо [math]y(x)[/math] рассматривается сеточная функция [math]\widehat{y}_{i}\approx y(x_{i})[/math], которая определяется только в точках сетки. 3. Выполняется подстановка аппроксимаций (6.16), (6.17) и [math]f(x_{i},y_{i})[/math] в дифференциальное уравнение: [math]\begin{gathered}\frac{y_{i+1}-y_{i}}{h_{i+1}}+O(h_{i+1})= f(x_{i},y_{i});\\ \frac{1}{H_{i}^{i+1}}\! \left[-\delta_{i+1}y_{i-1}+ \frac{\delta_{i+1}^2-1}{\delta_{i+1}}y_{i}+ \frac{y_{i+1}}{\delta_{i+1}}\right]+ O(h_{i}h_{i+1})= f(x_{i},y_{i}). \end{gathered}[/math] После отбрасывания остаточных слагаемых получаются разностные схемы: — явная схема Эйлера первого порядка (явный метод Эйлера): [math]\widehat{y}_{i+1}= \widehat{y}_{i}+ h_{i+1}\cdot f(x_{i},\widehat{y}_{i}),\quad i=\overline{0,n-1},\quad \widehat{y}_0=y_0;[/math] (6.19) – обобщенная на нерегулярный шаблон двухшаговая явная схема Эйлера второго порядка (где [math]\Delta \widehat{y}_{i}=\widehat{y}_{i+1}-\widehat{y}_{i}[/math]): [math]\widehat{y}_[i+1}= \widehat{y}_{i}-\delta_{i+1}^2 \Delta \widehat{y}_{i}+ H_{i}^{i+1} \delta_{i+1} f(x_{i},\widehat{y}_{i}),\quad i=\overline{1,n-1},[/math] (6.20) Обозначим схему (6.20) символами 2Y2B (суть этого обозначения поясняется далее). Анализ остаточного слагаемого [math]\left(\frac{h_{i}^2}{6}\delta_{i+1}M_{3,i}\right)[/math] обобщенной схемы Эйлера указывает на то, что при условном выборе шаблона, когда достигается соотношение [math]\delta_{i+1}\leqslant h_{i}[/math] (то есть [math]h_{i+1}\leqslant h_{i}^2[/math]) , порядок аппроксимации этой схемы повышается на единицу без увеличения количества точек шаблона. Этот прием повышения точности может быть использован при конструировании вложенных алгоритмов путем измельчения шага. При [math]\delta_{i+1}=1[/math] (регулярный шаблон) схема (6.20) соответствует методу Эйлера-Коши: [math]\widehat{y}_{i+1}= \widehat{y}_{i-1}+ 2h\cdot f(x_{i},\widehat{y}_{i}),\quad i=\overline{1,n-1}.[/math] (6.21) Для начала расчетов по формулам (6.20),(6.21) требуется иметь две "разгонные" точки [math]\widehat{y}_{0},\,\widehat{y}_{1}[/math]. Первая определяется известным начальным условием [math]\widehat{y}_0=y_0[/math], а вторая может быть найдена с помощью другого метода, например, по формуле (6.19): [math]\widehat{y}_1=y_0+h_1f(x_0,y_0)[/math]. Приведем геометрическую интерпретацию явного метода Эйлера (6.19), которая показана на рис. 6.5. Пусть известна точка [math](x_{i},\widehat{y}_{i})[/math] на искомой интегральной кривой [math]y=y(x)[/math]. Через эту точку проведем к кривой [math]y=y(x)[/math] касательную [math]L[/math], тангенс угла наклона которой равен [math]y'(x_{i})= f(x_{i}, \widehat{y}_{i})[/math]. Ее уравнение имеет вид [math]y=\widehat{y}_{i}+ y'(x_{i})(x-x_{i})= \widehat{y}_{i}+ f(x_{i}, \widehat{y}_{i})(x-x_{i}).[/math] Следующая точка [math]\widehat{y}_{i+1}[/math] приближенного решения получается при [math]x=x_{i+1}[/math] (в результате нахождения точки пересечения прямых [math]L[/math] и [math]x_{i+1}=\text{const}[/math]): [math]\widehat{y}_{i+1}= \widehat{y}_{i}+ h_{i+1}f(x_{i},\widehat{y}_{i})[/math]. Это соотношение совпадает с формулой (6.19). Схема (6.19) иногда называется схемой ломаных (рис. 6.6), так как очередное значение [math]\widehat{y}_{i+1}[/math] получается в результате проведения касательной к интегральной кривой в точке [math](x_{i},\widehat{y}_{i})[/math]. После выполнения каждого шага метода Эйлера приближенное решение переходит с одной интегральной кривой на другую, что обусловлено погрешностью численного решения. Большая погрешность схемы связана с тем, что порядок этой схемы равен единице, и поэтому ломаная, определяемая по (6.19), довольно значительно отклоняется от точного решения [math]y=y(x)[/math]. В связи с этим данная схема используется, как правило, в качестве вспомогательной при конструировании схем повышенного порядка. Для анализа устойчивости явного метода Эйлера применим формулу (6.19) при [math]h_{i+1}=h[/math] к решению задачи (6.12). В результате разностное уравнение будет иметь вид [math]\widehat{y}_{i+1}= \widehat{y}_{i}+ h\,\mu\,\widehat{y}_{i}[/math] или [math]\widehat{y}_{i+1}-(1+h\,\mu)\widehat{y}_{i}=0[/math]. Корень соответствующего характеристического уравнения [math]\lambda_1=1+h\,\mu[/math]. Условие (6.13) принимает форму [math]|1+h\,\mu|<1[/math]. Отсюда следует, что область устойчивости представляется внутренней частью круга единичного радиуса с центром в точке [math](-1;0)[/math] (рис. 6.7). Если [math]\mu[/math] действительная константа, то можно указать интервал устойчивости [math]-2<h\,\mu<0[/math], то есть [math]h\,\mu\in (-2;0)[/math], а при [math]\mu<0[/math] имеем [math]0<h<-\frac{2}{\mu}[/math]. Поэтому явный метод Эйлера является ограниченно устойчивым с критическим шагом [math]h_{\text{kr}}=-\frac{2}{\mu}[/math]. Дадим геометрическую интерпретацию метода Эйлера-Коши (6.21) применительно к решению задачи (6.9) (рис. 6.8). Пусть известны две точки [math](x_{i-1},\widehat{y}_{i-1}),\, (x_{i}, \widehat{y}_{i})[/math] на искомой интегральной кривой [math]y=y(x)[/math]. Проведем через точку [math](x_{i},\widehat{y}_{i})[/math] касательную [math]L[/math], тангенс угла наклона которой равен [math]y'(x_{i})= f(x_{i},\widehat{y}_{i})[/math]. После этого через точку [math](x_{i-1}, \widehat{y}_{i-1})[/math] проведем прямую [math]L_1[/math], параллельную [math]L[/math]. Ее уравнение имеет вид [math]y=\widehat{y}_{i-1}+ f(x_{i},\widehat{y}_{i})(x-x_{i-1})[/math]. Следующая точка [math]\widehat{y}_{i+1}[/math], приближенного решения получается при [math]x=x_{i+1}\colon[/math] [math]\widehat{y}_{i+1}= \widehat{y}_{i-1}+2h f(x_{i}, \widehat{y}_{i})[/math]. Это соотношение совпадает с формулой (6.21). Замечание. Все описанные здесь и далее методы легко переносятся на системы обыкновенных дифференциальных уравнений (6.7). Под [math]y[/math] и [math]f(x,y)[/math] следует подразумевать векторы [math]Y=(y_1,\ldots,y_n)^T[/math] и [math]F(x,Y)= \bigl(f_1(x,Y), \ldots, f_n(x,Y)\bigr)^T[/math]. Например, для системы двух дифференциальных уравнений [math]\begin{cases}y'=f(x,y,z),& y(x_0)=y_0,\\ z'=g(x,y,z),& z(x_0)=z_0 \end{cases}[/math] явный метод Эйлера (6.19) записывается в форме [math]\left\{\! \begin{aligned}&\widehat{y}_{i+1}= \widehat{y}_{i}+ h_{i+1}\cdot f(x,\widehat{y}_{i},\widehat{z}_{i}), &~ & \widehat{y}_0=y_0,\\ &\widehat{z}_{i+1}= \widehat{z}_{i}+ h_{i+1}\cdot g(x,\widehat{y}_{i},\widehat{z}_{i}), &~ & \widehat{z}_0=z_0. \end{aligned}\right.[/math] Пример 6.3. Исследовать устойчивость явного метода Эйлера на примере решения уравнения [math]Ty'+y=1,[/math] [math]y(0)=y_0[/math], где [math]T[/math] — действительное положительное число. Найти приближенное решение задачи Коши [math]0,\!1y'+y=1,~y(0)=0[/math] на отрезке [math][0;1][/math] явным методом Эйлера. РешениеТочное решение задачи получено в примере 6.1: [math]y(x)=1-e^{-10x}[/math], поскольку [math]T=0,\!1[/math]. Уравнение [math]Ty'+y=1[/math] перепишем в виде, аналогичном (6.12): [math]y'=-\frac{1}{T}y+\frac{1}{T}[/math], поэтому [math]\mu=-\frac{1}{T}[/math], а критический шаг равен [math]h_{\text{kr}}=-\frac{2}{\mu}=2T[/math]. Следовательно, при [math]h<2T[/math] метод устойчив, а при [math]h>2T[/math] неустойчив. На рис. 6.9 схематически изображены приближенные решения уравнения [math]Ty'+y=1,[/math] [math]y(0)=0[/math], соответствующие устойчивому процессу и хорошей точности численного решения (рис.6.9,а), устойчивому процессу и плохой точности (рис.6.9,б), а также неустойчивому процессу и плохой точности (рис. 6.9,в). Для уравнения [math]0,\!1y'+y=1[/math] постоянная времени [math]T=0,\!1[/math], поэтому [math]h_{\text{kr}}==0,\!2[/math]. Решим задачу для значений шага [math]h=0,\!05; 0,\!2; 0,\!5[/math]. При [math]h=0,\!05[/math] метод является устойчивым, при [math]h=0,\!2[/math] находится на фанице устойчивости, а при [math]h=0,\!5[/math] является неустойчивым. Переписав уравнение в виде [math]y'=10-10y[/math], получим из (6.19) формулу для нахождения приближенного решения: [math]\widehat{y}_{i+1}= \widehat{y}_{i}+ h\cdot (10-10\widehat{y}_{i})= 10h+ (1-10h)\widehat{y}_{i}.[/math] [math]\begin{array}{|c|c|c|c|c|} \multicolumn{5}{r}{\mathit{Table~6.1}}\\\hline x_{i}& \widehat{y}_{i}~ (h=0,\!05) & \widehat{y}_{i}~ (h=0,\!2) & \widehat{y}_{i}~ (h=0,\!5)& y(x_{i})\\\hline 0,\!00& 0,\!000000& 0,\!000& 0,\!000& 0,\!000000 \\\hline 0,\!05& 0,\!500000& & & 0,\!393469 \\\hline 0,\!01& 0,\!750000& & & 0,\!632121 \\\hline 0,\!15& 0,\!875000& & & 0,\!776870 \\\hline 0,\!20& 0,\!937500& 2,\!000& & 0,\!864665 \\\hline 0,\!25& 0,\!968750& & & 0,\!917915 \\\hline 0,\!30& 0,\!984375& & & 0,\!950213 \\\hline 0,\!35& 0,\!992187& & & 0,\!969803 \\\hline 0,\!40& 0,\!996094& 0,\!000& & 0,\!981684 \\\hline 0,\!45& 0,\!998047& & & 0,\!988891 \\\hline 0,\!50& 0,\!999023& & 5,\!000& 0,\!993262 \\\hline 0,\!55& 0,\!999511& & & 0,\!995913 \\\hline 0,\!60& 0,\!999755& 2,\!000& & 0,\!997521 \\\hline 0,\!65& 0,\!999877& & & 0,\!998497 \\\hline 0,\!70& 0,\!999938& & & 0,\!999088 \\\hline 0,\!75& 0,\!999969& & & 0,\!999446 \\\hline 0,\!80& 0,\!999984& 0,\!000& & 0,\!999665 \\\hline 0,\!85& 0,\!999992& & & 0,\!999797 \\\hline 0,\!90& 0,\!999996& & & 0,\!999877 \\\hline 0,\!95& 0,\!999998& & & 0,\!999925 \\\hline 1,\!00& 0,\!999999& 2,\!000& -15,\!000& 0,\!999955 \\\hline \max\varepsilon_{i}& 0,\!117879& 1,\!135335& 15,\!99995& \\\hline \end{array}[/math] Результаты расчетов, выполненных для различных значений [math]h[/math], приведены в табл. 6.1. При [math]h=0,\!05[/math] приближенное решение достаточно близко к точному, при [math]h=0,\!2[/math] "пила" не сходится и не расходится, а при [math]h=0,\!5[/math] "пила", аппроксимирующая точное решение, расходится (рис. 6.9,в). Замечание. Порядок точности метода, как правило, определяется порядком аппроксимации схем. Например, для схемы (6.19) первого порядка, записанной в виде, соответствующем исходному уравнению (6.9), получаем (первое слагаемое аппроксимирует производную): [math]\frac{y(x_{i+1})-y(x_{i})}{h_{i+1}}-f(x_{i},y_{i})= \frac{y_{i}+ h_{i+1}y'_{i}+ O(h_{i+1}^2)-y_{i}}{h_{i+1}}-f_{i}= y'_{i}-f_{i}+ O(h_{i+1})=0,[/math] то есть [math]y'_{i}=f_{i}+O(h_{i+1})[/math]. Здесь были использованы формула Тейлора и свойства, приведенные ранее. Для двухшаговой схемы (6.21) аналогично получается второй порядок аппроксимации. Действительно, воспользуемся частным случаем формулы Тейлора (см. (В. 18)) при [math]x_{0}=x_{i},\quad x_{i+1}=x_{i}+h,\quad x_{i-1}=x_{i}-h_{i},\quad f(x)= y(x)[/math] имеем [math]y_{i+1}= y_{i}+ hy'_{i}+ \frac{h^2}{2}y''_{i}+ O(h^3),\qquad y_{i-1}= y_{i}-hy'_{i}+ \frac{h^2}{2}y''_{i}+ O(h^3).[/math] Подставляя полученные соотношения в схему (6.21), записанную в виде [math]\frac{y_{i+1}-y_{i-1}}{2h}-f(x_{i},y_{i})=0[/math], соответствующем дифференциальному уравнению (6.9), с учетом свойств, приведенных ранее, имеем [math]\frac{y_{i}+ hy'_{i}+ \dfrac{h^2}{2}y''_{i}+ O(h^3)-y_{i}+ hy'_{i}-\dfrac{h^2}{2}y''_{i}-O(h^3)}{2h}-f_{i}= y'_{i}+ O(h^2)-f_{i}=0[/math], то есть [math]y'_{i}=f_{i}+O(h^2)[/math]. Также легко проверить, что и схема (6.20) на нерегулярном шаблоне имеет второй порядок при безусловной аппроксимации. Если же [math]h_{i+1}\leqslant h_{i}^2[/math] (условная аппроксимация), то этот порядок равен трем. Интегрально-интерполяционный принципВ данном подходе явные или неявные методы соответствующего порядка получаются путем интегрирования уравнения (6.9) вдоль решения [math]y(x)[/math]. Так, на двухточечном шаблоне [math]\int\limits_{y_{i}}^{y_{i+1}}dy= \int\limits_{x_{i}}^{x_{i+1}}f(x,y(x))\,dx\,.[/math] (6.22) Функция [math]f(x,y(x))[/math] в (6.22) заменяется интерполяционным многочленом соответствующей степени [math]p[/math] ([math]p=0[/math] для метода первого порядка, [math]p=1[/math] для метода второго порядка, [math]p=2[/math] для метода третьего порядка и т.д.). Подчеркнем, что замена [math]f(x,y)[/math] интерполяционным многочленом в (6.22) становится возможной благодаря тому, что [math]f(x,y)[/math] рассматривается не на всей плоскости [math](x,y)[/math], а только вдоль интегральной кривой [math]y(x)[/math]. Для вычисления интеграла в правой части могут быть использованы известные квадратурные формулы прямоугольников, трапеций, парабол (Симпсона) и др. Примем [math]x_{i+1}[/math] в качестве узловой точки многочлена [math]L_{0}(x)[/math], то есть [math]L_0(x)= f(x_{i+1},\widehat{y}_{i+1})[/math]. Тогда из (6.22) следует формула неявного метода Эйлера: [math]\widehat{y}_{i+1}= \widehat{y}_{i}+ h_{i+1} f(x_{i+1}, \widehat{y}_{i+1})\equiv \Phi(x_{i},x_{i+1},\widehat{y}_{i+1}),\quad i=\overline{0,n-1}.[/math] (6.23) Принимая для интерполяционного многочлена [math]L_1(x)[/math] в качестве узловых точек [math]x_{i}[/math] и [math]x_{i+1}[/math] и используя квадратурную формулу трапеций (5.46), получаем неявную одношаговую схему второго порядка (метод трапеций): [math]\widehat{y}_{i+1}= \widehat{y}_{i}+ \frac{h_{i+1}}{2}\bigl[f_{i}+ f(x_{i+1}, \widehat{y}_{i+1})\bigr]\equiv \Phi(x_{i}, x_{i+1}, \widehat{y}_{i+1}),\quad i=\overline{0,n-1}.[/math] (6.24) где [math]f_{i}= f(x_{i},\widehat{y}_{i+1})[/math]. Формулы (6.23),(6.24) имеют вид (6.11) при [math]k=1[/math]. Подчеркнем, что свойство неявности схем (6.23),(6.24) обусловлено наличием искомой величины [math]\widehat{y}_{i+1}[/math] в левой и правой частях в общем случае нелинейных уравнений. При реализации алгоритма решения задачи Коши неизвестное значение [math]\widehat{y}_{i+1}[/math] вычисляется одним из методов решения нелинейных уравнений. Применение метода Ньютона связано с записью уравнения в форме [math]\widehat{y}_{i+1}-\Phi(x_{i},x_{i+1}, \widehat{y}_{i+1})\equiv F(\widehat{y}_{i+1})=0[/math] (6.25) и с дифференцированием функции [math]F(\widehat{y}_{i+1})[/math], что увеличивает время расчетов из-за возможной сложности производных или увеличивает погрешность в силу неточного задания правых частей. Как правило, используется метод простых итераций: [math]\widehat{y}_{i+1}^{\,(k+1)}= \Phi \bigl(x_{i}, x_{i+1}, \widehat{y}_{i+1}^{\,(k)}\bigr),\quad k=0,1,2,\ldots[/math] (6.26) Таким образом, уравнение (6.24) решается с помощью формулы (рекуррентного соотношения): [math]\widehat{y}_{i+1}^{\,(k+1)}= \widehat{y}_{i}+ \frac{h_{i+1}}{2} \bigl[f_{i}+ f(x_{i+1}, \widehat{y}_{i+1}^{\,(k)})\bigr],\quad k=0,1,2,\ldots[/math] (6.27) При применении методов Ньютона и простых итераций вначале задается или находится нулевое приближение решения по формуле [math]y_{i+1}^{\,(0)}= \widehat{y}_{i}[/math], (так называемый "постоянный" прогноз) или явным методом Эйлера (6.19): [math]\widehat{y}_{i+1}^{\,(0)}= \widehat{y}_{i}+ h_{i+1} f(x_{i}, \widehat{y}_{i})[/math]. Итерации завершаются при выполнении условия окончания [math]\bigl|\widehat{y}_{i+1}^{\,(k+1)}-\widehat{y}_{i+1}^{\,(k)}\bigr| \leqslant \varepsilon[/math], где [math]\varepsilon[/math] — малое положительное число. Для сокращения описания и удобства ссылок получаемые здесь схемы обозначаются буквенно-цифровыми символами. Первая цифра в нем обозначает "шаговость" схемы, следующая за ней буква "Я" или буквы "НЯ" — явный или неявный характер, затем цифра указывает порядок точности схемы и за ней буква — возможную модификацию. Таким образом, схема (6.24) обозначается 1НЯ2, т.е. она является одношаговой, неявной и имеет второй порядок точности. Последняя буква указывается только для тех схем, которые имеют модификации. Приведем геометрическую интерпретацию неявного метода Эйлера (6.23) (рис.6.10). Пусть известна точка [math](x_{i}, \widehat{y}_{i})[/math] на искомой интегральной кривой [math]y=y(x)[/math]. Через эту точку проведем прямую, тангенс угла наклона которой равен [math]y'(x_{i+1})= f(x_{i+1}, \widehat{y}_{i+1})[/math], т.е. следующая точка [math]\widehat{y}_{i+1}[/math] приближенного решения берется на той кривой из семейства интегральных кривых, касательная к которой в точке [math](x_{i+1},\widehat{y}_{i+1})[/math] проходит через точку [math](x_{i},\widehat{y}_{i})[/math]. Поэтому [math]\frac{\widehat{y}_{i+1}-\widehat{y}_{i}}{h_{i+1}}= f(x_{i+1}, \widehat{y}_{i+1})[/math]. Полученное соотношение совпадает с формулой (6.23). Исследуем устойчивость неявного метода Эйлера и метода трапеций. Применим формулу (6.23) при [math]h_{i+1}=h[/math] к решению задачи (6.12). В результате получим разностное уравнение [math]\widehat{y}_{i+1}= \widehat{y}_{i}+h\cdot\mu\cdot \widehat{y}_{i+1}[/math] или [math]\widehat{y}_{i+1}(1-h\cdot\mu)-\widehat{y}_{i}=0[/math]. Корень характеристического уравнения [math]\lambda_1= \frac{1}{1-h\cdot\mu}[/math]. Условие устойчивости (6.13) имеет вид [math]\left|\frac{1}{1-h\cdot\mu}\right|<1[/math] или [math]|1-h\cdot\mu|>1[/math]. Область устойчивости представляется внешней частью круга единичного радиуса с центром в точке [math](1;0)[/math]. Из геометрической интерпретации (рис. 6.11) следует, что неявный метод Эйлера обладает свойством A-устойчивости (см. рис.6.4,а), поскольку выделенная на рис. 6.11 область включает левую полуплоскость. Применим формулу (6.24) метода трапеций при [math]h_{i+1}=h[/math] к решению задачи (6.12). Получим разностное уравнение [math]\widehat{y}_{i+1}= \widehat{y}_{i}+ \frac{h}{2} (\mu \widehat{y}_{i}+ \mu \widehat{y}_{i+1})[/math] или [math]\widehat{y}_{i+1}\! \left(1-\frac{h\mu}{2}\right)-\left(1+ \frac{h\mu}{2}\right)\! \widehat{y}_{i}=0[/math]. Корень характеристического уравнения [math]\lambda_1= \left(1+ \frac{h\mu}{2}\right)\,\colon \left(1-\frac{h\mu}{2}\right)[/math]. Условие устойчивости (6.13) принимает вид [math]\left|\left(1+ \frac{h\mu}{2}\right)\,\colon \left(1-\frac{h\mu}{2}\right)\right|<1[/math] или [math]|2+m\mu|< |2-h\mu|[/math]. Полученное неравенство выполняется в области, изображенной на рис. 6.4,а, т.е. метод трапеций является A-устойчивым. Продолжим рассмотрение интегрально-интерполяционного метода. Использование квадратурной формулы прямоугольников дает модифицированный метод Эйлера второго порядка: [math]\widehat{y}_{i+1\!\not{\phantom{|}}\,\,2}= \widehat{y}_{i}+ \frac{h_{i+1}}{2} f(x_{i}, \widehat{y}_{i}),\quad i=\overline{0,n-1},[/math] [math]\widehat{y}_{i+1}= \widehat{y}_{i}+ h_{i+1} f\! \left(x_{i}+ \frac{h_{i+1}}{2}, \widehat{y}_{i+1\!\not{\phantom{|}}\,\,2}\right)\!,\quad i=\overline{0,n-1}.[/math] (6.28) Дадим геометрическую интерпретацию модифицированного метода Эйлера второго порядка (рис. 6.12). Пусть известна точка [math](x_{i},\widehat{y}_{i})[/math] на искомой интегральной кривой [math]y=y(x)[/math]. Через эту точку проведем касательную [math]L[/math], тангенс угла наклона которой равен [math]y'(x_{i})= f(x_{i},\widehat{y}_{i})[/math]. При [math]x=x_{i}+ \frac{h_{i+1}}{2}[/math] получаем точку [math]\widehat{y}_{i+1\!\not{\phantom{|}}\,\,2}[/math]. Далее проведем прямую [math]L_1[/math] с тангенсом угла наклона [math]f\! \left(x_{i}+ \frac{h_{i+1}}{2}, \widehat{y}_{i+1\!\not{\phantom{|}}\,\,2}\right)[/math], а затем прямую [math]L_2[/math], параллельную [math]L_1[/math] и проходящую через точку [math](x_{i}, \widehat{y}_{i})[/math]. Ее уравнение имеет вид [math]y=\widehat{y}_{i}+ f\! \left(x_{i}+ \frac{h_{i+1}}{2}, \widehat{y}_{i+1\!\not{\phantom{|}}\,\,2}\right)\!\cdot (x-x_{i})[/math]. Следующая точка [math]\widehat{y}_{i+1}[/math] приближенного решения находится при [math]x=x_{i+1}[/math] (это соотношение совпадает с формулой (6.28)): [math]\widehat{y}_{i+1}= \widehat{y}_{i}+ h_{i+1} f\! \left(x_{i}+ \frac{h_{i+1}}{2}, \widehat{y}_{i+1\!\not{\phantom{|}}\,\,2}\right)\!.[/math] Исследуем устойчивость модифицированного метода Эйлера второго порядка. Применим формулы (6.28) при [math]h_{i+1}=h[/math] к решению задачи (6.12). В результате получим разностное уравнение: [math]\widehat{y}_{i+1}= \widehat{y}_{i}+ h\mu \left(\widehat{y}_{i}+ \frac{h}{2}\mu \widehat{y}_{i}\right)[/math] или [math]\widehat{y}_{i+1}-\widehat{y}_{i} \left(1+ h\mu+ \frac{(h\mu)^2}{2}\right)=0[/math]. Корень соответствующего характеристического уравнения [math]\lambda_1= 1+h\mu+ \frac{(h\mu)^2}{2}[/math]. Условие устойчивости (6.13) имеет вид [math]\left|1+h\mu+ \frac{(h\mu)^2}{2}\right|<1[/math]. Если [math]\mu[/math] — действительное число, то, раскрывая модуль, получаем [math]-1<1+h\mu+ \frac{(h\mu)^2}{2}<1[/math]. Неравенство [math]2+h\mu+ \frac{(h\mu)^2}{2}>0[/math] выполняется для любых [math]h\mu[/math], а неравенство [math]h\mu+ \frac{(h\mu)^2}{2}<0[/math], то есть [math]h\mu\cdot (2+h\mu)<0[/math], справедливо при [math]-2<h\mu<0[/math]. Таким образом, интервал устойчивости [math]h\mu\in (-2;0)[/math] модифицированного метода Эйлера совпадает с интервалом устойчивости явного метода Эйлера (6.19). Для использования регулярного трехточечного шаблона [math](x_{i-1},x_{i},x_{i+1})[/math] заменим (6.22) аналогичным выражением: [math]\textstyle{\int\limits_{y_{i-1}}^{y_{i+1}} dy= \int\limits_{x_{i-1}}^{x_{i+1}} f(x,y(x))dx}[/math]. Если функцию [math]f(x,y(x))[/math], входящую в правую часть, заменить многочленом Лагранжа второй степени и применить формулу (5.47), то получается неявная схема парабол (Симпсона) четвертого порядка (2НЯ4): [math]\widehat{y}_{i+1}= \widehat{y}_{i-1}+ \frac{h}{3} \bigl[f_{i-1}+ 4f_{i}+ f(x_{i+1}, \widehat{y}_{i+1})\bigr].[/math] (6.29) Отметим, что схема парабол имеет четвертый, а не третий порядок точности из-за регулярности шаблона. В соответствии с классификацией, схема (6.29) относится к классу неявных схем Адамса. Для получения других схем Адамса следует переписать (6.22) в форме [math]\int\limits_{y_{i-1}}^{y_{i+1}} dy= \int\limits_{x_{i-1}}^{x_{i+1}} f(x,y(x))\,dx= \int\limits_{x_{i}}^{x_{i+1}} y'(x)\,dx[/math] (6.30) и заменить [math]y'(x)[/math] соответствующим интерполяционным полиномом, построенным на регулярном шаблоне с [math]h=\text{const}[/math]. В силу формулы (4.36) второго интерполяционного многочлена Ньютона получаем (где [math]\widehat{q}= \frac{x-x_{i}}{h}[/math] — фаза интерполяции) [math]y'(\widehat{q})= y'_{i}+ \widehat{q} \Delta y'_{i-1}+ \frac{\Delta^2 y'_{i-2}}{2!} \widehat{q} (\widehat{q}+1)+ \frac{\Delta^3 y'_{i-3}}{3!} \widehat{q} (\widehat{q}+1)(\widehat{q}+2)+\ldots,[/math] или [math]y'(\widehat{q})= y'_{i}+ \widehat{q} \Delta y'_{i-1}+ \frac{\widehat{q}^2+ \widehat{q}}{2} \Delta^2 y'_{i-2}+ \frac{\widehat{q}^3+ 3\widehat{q}^2+ 2\widehat{q}}{6} \Delta^3 y'_{i-3}+\ldots[/math] Подставляя полученное выражение в правую часть (6.30) и учитывая, что [math]dx= h\cdot d\widehat{q}[/math], вычисляем [math]\textstyle{\int\limits_{0}^{1} y'(\widehat{q})\cdot h\, d\widehat{q}}[/math], поскольку при [math]x=x_{i}[/math] и [math]x=x_{i+1}[/math] справедливо [math]\widehat{q}=0[/math] и [math]\widehat{q}=1[/math] соответственно. В результате имеем [math]\widehat{y}_{i+1}-\widehat{y}_{i}= h\cdot y'_{i}+ \frac{1}{2}h\cdot \Delta y'_{i-1}+ \frac{5}{12}h\cdot \Delta^2 y'_{i-2}+ \frac{3}{8}h\cdot \Delta^3 y'_{i-3}+\ldots[/math] или с учетом (6.9) общую формулу семейства явных схем Адамса: [math]\widehat{y}_{i+1}= \widehat{y}_{i}+ \eta_{i}+ \frac{1}{2} \Delta \eta_{i-1}+ \frac{5}{12} \Delta^2 \eta_{i-2}+ \frac{3}{8} \Delta^3 \eta_{i-3}+ \frac{251}{720} \Delta^4 \eta_{i-4}+\ldots,[/math] (6.31) где [math]\eta_{i}= h\cdot f(x_{i}, \widehat{y}_{i})= h\cdot f_{i}[/math], а конечные разности вычисляются с помощью табл. 6.2. Подчеркнем, что здесь интерполяционный многочлен записан на шаблоне, не включающем точку [math]x_{i+1}[/math]. Если взять в правой части (6.31) два слагаемых, то получится явный метод Эйлера (6.19), если три слагаемых, то явный метод, описываемый формулой [math]\widehat{y}_{i+1}= \widehat{y}_{i}+ h\cdot f_{i}+ \frac{1}{2} (h\cdot f_{i}-h\cdot f_{i-1})= \widehat{y}_{i}+ \frac{h}{2} (3f_{i}-f_{i-1}).[/math] а если четыре слагаемых, то явный метод вида [math]\begin{aligned}\widehat{y}_{i+1}&= \widehat{y}_{i}+ \frac{h}{2} (3f_{i}-f_{i-1})+ \frac{5}{12} \bigl[(h\cdot f_{i}-h\cdot f_{i-1})-(h\cdot f_{i-1}-h\cdot f_{i-2})\bigr]=\\ &=\widehat{y}_{i}+ \frac{h}{12} \bigl(23f_{i}-16f_{i-1}+ 5f_{i-2}\bigr) \end{aligned}[/math] В результате имеем схемы Адамса-Бэшфорта – второго порядка (2Я2В): [math]\widehat{y}_{i+1}= \widehat{y}_{i}+ \frac{h}{2} (3f_{i}-f_{i-1}),\quad i=\overline{1,n-1};[/math] (6.32) – третьего порядка (ЗЯЗА): [math]\widehat{y}_{i+1}= \widehat{y}_{i}+ \frac{h}{12} \bigl(23f_{i}-15f_{i-1}+ 5f_{i-2}\bigr),\quad i=\overline{2,n-1};[/math] (6.33) – четвертого порядка (4Я4): [math]\widehat{y}_{i+1}= \widehat{y}_{i}+ \frac{h}{24} \bigl(55f_{i}-59f_{i-1}+ 37f_{i-2}-9f_{i-3}\bigr),\quad i=\overline{3,n-1}.[/math] (6.34) Для начала расчетов по формуле (6.32) требуются две "разгонные" точки: [math]\widehat{y}_{0},\, \widehat{y}_{1}[/math], по формуле (6.33) — три "разгонные" точки: [math]\widehat{y}_{0},\, \widehat{y}_{1},\, \widehat{y}_{2}[/math], по формуле (6.34) — четыре "разгонные" точки: [math]\widehat{y}_{0},\, \widehat{y}_{1},\, \widehat{y}_{2},\, \widehat{y}_{3}[/math]. Их необходимо вычислить с порядком точности не меньше порядка точности схемы. Далее схема (6.32) обобщается на нерегулярный шаблон. Если для построения интерполяционного многочлена использовать шаблоны, включающие точку [math]x_{i+1}[/math], то получится следующая формула семейства неявных схем Адамса: [math]\widehat{y}_{i+1}= \widehat{y}_{i}+\eta_{i+1}-\frac{1}{2} \Delta \eta_{i}-\frac{1}{12} \Delta^2\eta_{i-1}-\frac{1}{24} \Delta^3\eta_{i-2}-\frac{19}{720} \Delta^2 \eta_{i-3}+\ldots,[/math] (6.35) где [math]\eta_{i+1}= h\cdot f(x_{i+1}, \widehat{y}_{i+1})= h\cdot f_{i+1}[/math], а конечные разности вычисляются с помощью табл. 6.3. Если взять в правой части (6.35) два слагаемых, то получится неявный метод Эйлера (6.23), если три слагаемых, — метод трапеций (6.24): [math]\widehat{y}_{i+1}= \widehat{y}_{i}+ h\cdot f_{i+1}-\frac{1}{2} (h\cdot f_{i+1}-h\cdot f_{i})= \widehat{y}_{i}+ \frac{h}{2} \bigl[f_{i}+ f(x_{i+1}, \widehat{y}_{i+1})\bigr],[/math] если четыре слагаемых, — неявный метод, описываемый соотношением (и т.д.) [math]\begin{aligned}\widehat{y}_{i+1}&= \widehat{y}_{i}+ \frac{h}{2} (f_{i+1}+f_{i+1})-\frac{1}{12} \bigl[(h\cdot f_{i+1}-h\cdot f_{i})-(h\cdot f_{i}-h\cdot f_{i-1})\bigr]=\\ &=\widehat{y}_{i}+ \frac{h}{12} \bigl[-f_{i-1}+ 8f_{i}+ 5f(x_{i+1}, \widehat{y}_{i+1})\bigr]. \end{aligned}[/math] В результате получаются неявные схемы Адамса-Мултона: – первого порядка (6.23); – второго порядка (6.24); – третьего порядка (2НЯЗ): [math]\widehat{y}_{i+1}= \widehat{y}_{i}+ \frac{h}{12} \bigl[-f_{i-1}+ 8f_{i}+ 5f(x_{i+1}, \widehat{y}_{i+1})\bigr],\quad i=\overline{1,n-1};[/math] (6.36) – четвертого порядка (ЗНЯ4): [math]\widehat{y}_{i+1}= \widehat{y}_{i}+ \frac{h}{24} \bigl[f_{i-2}-5f_{i-1}+ 19f_{i}+ 9f(x_{i+1}, \widehat{y}_{i+1})\bigr],\quad i=\overline{2,n-1}.[/math] (6.37) Для расчетов по формуле (6.36) требуются две "разгонные" точки: [math]\widehat{y}_{0},\, \widehat{y}_{1}[/math], по формуле (6.37) — три "разгонные" точки: [math]\widehat{y}_{0},\, \widehat{y}_{1},\, \widehat{y}_{2}[/math]. Их необходимо вычислить с порядком точности не меньше порядка точности схемы. Чтобы найти искомое значение [math]\widehat{y}_{i+1}[/math] с помощью формул (6.36),(6.37), так же как в неявном методе Эйлера и методе трапеций, требуется решить в общем случае нелинейное уравнение. Замечание. Аналогичные методы могут быть получены из интегрального уравнения вида [math]\textstyle{\int\limits_{y_{i-k}}^{y_{i+1}} dy= \int\limits_{x_{i-k}}^{x_{i+1}} f(x,y(x))\,dx,~ k \geqslant 2}[/math]. ▼ Примеры 6.4-6.5
Принцип согласования с разложением по формуле ТейлораРассмотрим задачу Коши (6.9) в общем случае на неравномерной сетке [math]\Omega_{n}[/math]. Для получения формулы нахождения последующего значения [math]\widehat{y}_{i+1}[/math] на шаблоне [math](x_{i}, x_{i+1})[/math] разложим решение [math]y(x)[/math] при [math]x=x_{i+1}[/math] относительно точки [math]x_{i}[/math] по формуле Тейлора на промежутке [math]x_{i} \leqslant x \leqslant x_{i+1}\colon[/math] [math]y_{i+1}= y_{i}+ h_{i+1}y'_{i}+ \frac{1}{2}h_{i+1}^2y''_{i}+ \ldots+ \frac{1}{n!} h_{i+1}^{p} y_{i}^{(p)}+ O(h_{i+1}^{p+1}).[/math] (6.38) Входящие в правую часть производные можно определить, дифференцируя уравнение (6.9) требуемое число раз: [math]y'=f(x,y),\quad y''= \frac{d}{dx}f(x,y)= f_x+f_y\cdot y'= f_x+f\cdot f_y,\quad \ldots[/math] Если функция [math]f(x,y)[/math] имеет p-е непрерывные производные по обоим аргументам, то в разложении (6.38) можно учесть слагаемые до [math]O(h_{i+1}^{p+1})[/math]. Чем больше членов разложения используется для вычисления, тем точнее будет приближение и тем выше порядок схемы. Однако рассчитывать производные путем дифференцирования правой части уравнения (6.9) невыгодно. Поэтому Рунге и Кутта независимо друг от друга предложили более удобную и рациональную форму представления численного решения, при использовании которой дифференцировать функцию [math]f(x,y)[/math] не требуется. В общем случае методы Рунге-Кутты имеют следующую структуру (где [math][/math]): [math]\widehat{y}_{i+1}= \widehat{y}_{i}+ h_{i+1}\cdot \bigl[b_1K_{1,i}+ b_2K_{2,i}+ \ldots+ b_sK_{s,i}\bigr],\quad \widehat{y}_{0}= y_0,\quad i=\overline{0,n-1},[/math] где [math]\begin{aligned}& K_{1,i}= f(x_{i},\widehat{y}_{i});\\ & K_{2,i}= f(x_{i}+ c_2h_{i+1},\, \widehat{y}_{i}+ h_{i+1}a_{2,1}K_{1,i});\\ & K_{3,i}= f(x_{i}+ c_3h_{i+1},\, \widehat{y}_{i}+ h_{i+1}~(a_{3,1}K_{1,i}+ a_{3,2}K_{2,i});\\ &\quad\vdots\\ & K_{s,i}= f(x_{i}+ c_{s}h_{i+1},\, \widehat{y}_{i}+ h_{i+1}~ (a_{s,1}K_{1,i}+ a_{s,2} K_{2,i}+ \ldots+ a_{s,s-1} K_{s-1,i}), \end{aligned}[/math] (6.39) где [math]s[/math] — число стадий (этапов), [math]K_{s,i}[/math] — значения коэффициентов схемы Рунге-Кутты, вычисленные на основе правой части уравнения (6.9), [math]c_{j},~ j= \overline{s,s};[/math] [math]a_{\ell,m},~ \ell= \overline{2,s};[/math] [math]m=\overline{1,s-1};[/math] [math]b_k,~ k=\overline{1,s}[/math]. Первый индекс в обозначениях коэффициентов является порядковым номером, а второй соответствует индексу точки [math]x_{i}[/math] — началу отрезка [math][x_{i},x_{i+1}][/math], на котором производится расчет. В некоторых методах кроме вычисления приближенного решения [math]\widehat{y}_{i+1}[/math] определяется еще дополнительное значение [math]\widetilde{y}_{i+1}[/math] по формуле [math]\widetilde{y}_{i+1}= \widehat{y}_{i}+ h_{i+1}\cdot \bigl[\widetilde{b}_{1} K_{1,i}+ \widetilde{b}_2 K_{2,i}+ \ldots+ \widetilde{b}_{s} K_{s,i}\bigr],[/math] (6.40) порядок которого, как правило, на единицу больше или меньше обеспечиваемого выражением для [math]\widehat{y}_{i+1}[/math]. Величина [math]|\widehat{y}_{i+1}-\widetilde{y}_{i+1}|[/math] служит для учета погрешности и управления величиной шага. Коэффициенты схемы подбираются так, чтобы выражение для [math]\widehat{y}_{i+1}[/math] согласовывалось с формулой Тейлора (6.38) вплоть до членов степени [math]h_{i+1}^{p}[/math], где степень [math]p[/math] меняется в зависимости от метода и соответствует порядку схемы. Наибольшее распространение в вычислительной практике нашел метод Рунге-Кутты четвертого порядка: [math]\widehat{y}_{i+1}= \widehat{y}_{i}+ \frac{h_{i+1}}{6} (K_{1,i}+ 2K_{2,i}+ 2K_{3,i}+ K_{4,i}),\quad \widehat{y}_{0}= y_0,\quad i=\overline{0,n-1},[/math] (6.41) где [math]K_{1,i}= f_{i}= f(x_{i},\widehat{y}_{i}),~~ K_{2,i}= f\! \left(x_{i}+ \frac{h_{i+1}}{2},\, \widehat{y}_{i}+ \frac{h_{i+1}}{2}K_{1,i}\right)\!,[/math] [math]K_{3,i}= f\! \left(x_{i}+ \frac{h_{i+1}}{2},\, \widehat{y}_{i}+ \frac{h_{i+1}}{2} K_{2,i}\right)\!,~~ K_{4,i}= f \bigl(x_{i}+ h_{i+1},\, \widehat{y}_{i}+ h_[i+1}\cdot K_{3,i}\bigr).[/math] Схема (6.41) является четырехчленной, первый коэффициент [math]K_{1,i}[/math] относится к точке [math]x_{i}[/math], второй и третий — к средней точке [math]x_{i}+ \frac{h_{i+1}}{2}[/math] четвертый — к точке [math]x_{i+1}[/math]. Для этой схемы выбираются следующие параметры, входящие в (6.39): [math]s=4,\quad c_2=c_3=\frac{1}{2};\quad c_4=1;\quad a_{2,1}= \frac{1}{2};\quad a_{3,1}= a_{4,1}= a_{4,2}=0;\quad a_{3,2}=\frac{1}{2};\quad a_{4,3}=1.[/math] Приведем геометрическую интерпретацию метода Рунге-Кутты четвертого порядка (рис. 6.13). Пусть известна точка [math](x_{i},\widehat{y}_{i})[/math] на искомой интегральной кривой [math]y=y(x)[/math]. Прямая [math]AE[/math], проходящая через точку [math](x_{i}, \widehat{y}_{i})[/math] и определяющая величину [math]\widehat{y}_{i+1}[/math], имеет угловой коэффициент. Этот коэффициент рассчитывается по угловым коэффициентам касательных в точках [math]A,B,C,D[/math], проведенным к проходящим через эти точки интегральным кривым (штриховые линии). При этом коэффициент [math]K_{1,i}[/math] рассчитывается в точке [math]A[/math], коэффициент [math]K_{2,i}[/math] — в точке [math]C[/math], коэффициент [math]K_{3,i}[/math] — в точке [math]B[/math], а коэффициент [math]K_{4,i}[/math] — в точке [math]D[/math]. Замечания 1. Схемы Рунге-Кутты легко обобщаются на системы дифференциальных уравнений. Так, применительно к системе двух уравнений [math]\begin{cases}y'(x)= f(x,y,z),& y(x_0)=y_0,\\ z'(x)= g(x,y,z),& z(x_0)=z_0 \end{cases}[/math] схема Рунге-Кутты четвертого порядка имеет вид [math](i=\overline{0,n-1})[/math] [math]\begin{aligned} &\widehat{y}_{i+1}= \widehat{y}_{i}+ \frac{h_{i+1}}{6} (K_{1,i}+ 2K_{2,i}+ 2K_{3,i}+ K_{4,i}), &\quad & \widehat{y}_{0}=y_0,\\ &\widehat{z}_{i+1}= \widehat{z}_{i}+ \frac{h_{i+1}}{6} (q_{1,i}+ 2q_{2,i}+ 2q_{3,i}+ q_{4,i}), &\quad & \widehat{z}_0=z_0, \end{aligned}[/math] (6.42) где [math]\begin{aligned}& K_{1,i}= f(x_{i}, \widehat{y}_{i}, \widehat{z}_{i}); &\quad & K_{2,i}= f\! \left(x_{i}+ \frac{h_{i+1}}{2},\, \widehat{y}_{i}+ \frac{h_{i+1}}{2} K_{1,i},\, \widehat{z}_{i}+ \frac{h_{i+1}}{2} q_{1,i}\right)\!;\\[2pt] & qK_{1,i}= g(x_{i}, \widehat{y}_{i}, \widehat{z}_{i}); &\quad & q_{2,i}= g\! \left(x_{i}+ \frac{h_{i+1}}{2},\, \widehat{y}_{i}+ \frac{h_{i+1}}{2} K_{1,i},\, \widehat{z}_{i}+ \frac{h_{i+1}}{2} q_{1,i}\right)\!;\end{aligned}[/math] [math]\begin{aligned}& K_{3,i}= f\! \left(x_{i}+ \frac{h_{i+1}}{2},\, \widehat{y}_{i}+ \frac{h_{i+1}}{2} K_{2,i},\, \widehat{z}_{i}+ \frac{h_{i+1}}{2} q_{2,i}\right)\!;\\[2pt] & q_{3,i}= g\! \left(x_{i}+ \frac{h_{i+1}}{2},\, \widehat{y}_{i}+ \frac{h_{i+1}}{2} K_{2,i},\, \widehat{z}_{i}+ \frac{h_{i+1}}{2} q_{2,i}\right)\!;\\[2pt] & K_{4,i}= f \bigl(x_{i}+h_{i+1},\, \widehat{y}_{i}+ h_{i+1}K_{3,i},\, \widehat{z}_{i}+ h_{i+1} q_{3,i}\bigr);\\[2pt] & q_{4,i}= g \bigl(x_{i}+h_{i+1},\, \widehat{y}_{i}+ h_{i+1}K_{3,i},\, \widehat{z}_{i}+ h_{i+1} q_{3,i}\bigr). \end{aligned}[/math] 2. Соотношения модифицированного метода Эйлера (6.28) можно записать в форме схемы Рунге-Кутты: [math]\widehat{y}_{i+1}= \widehat{y}_{i}+ h_{i+1}K_{2,i},\quad i=\overline{0,n-1},\quad \widehat{y}_{0}=y_0,[/math] (6.43) [math]K_{1,i}= f(x,\widehat{y}_{i}),\quad K_{2,i}= f\! \left(x_{i}+ \frac{h_{i+1}}{2},\, \widehat{y}_{i}+ \frac{h_{i+1}}{2}K_{1,i}\right)\!.[/math] Покажем, что эта схема второго порядка. Разложение (6.38) с учетом равенства [math]y''= f_x+ f_{i}\cdot f_{y}[/math], где [math]f_{i}= f(x_{i},y_{i})[/math], частные производные [math]f_x, f_y[/math] и вторая производная [math]y''[/math] вычисляются в точке [math](x_{i}, y_{i})[/math], принимает вид [math]y_{i+1}= y_{i}+ h_{i+1}f_{i}+ \frac{h_{i+1}^2}{2} (f_x+ f_{i}\cdot f_y)+ O(h_{i+1}^3).[/math] (6.44) Разложение функции [math]f(x,y)[/math] относительно точки [math](x_{i},y_{i})[/math] можно записать следующим образом: [math]f(x,y)= f(x_{i},y_{i})+ (x-x_{i})f_x+ (y-y_{i})f_y+\ldots[/math] (6.45) где частные производные [math]f_x,f_y[/math] вычисляются в точке [math](x_{i},y_{i})[/math]. Тогда при [math]x=x_{i}+\frac{h_{i+1}}{2},~ y=y_{i}+\frac{h_{i+1}}{2}f_{i}[/math] получаем [math]f\! \left(x_{i}+\frac{h_{i+1}}{2},\, y_{i}+\frac{h_{i+1}}{2}f_{i}\right)= f_{i}+ \frac{h_{i+1}}{2}f_x+ \frac{h_{i+1}}{2}f_{i}\cdot f_y+ O(h_{i+1}^2).[/math] Следовательно, формулу (6.43) можно переписать в виде [math]\widehat{y}_{i+1}= \widehat{y}_{i}+ h_{i+1} f\! \left(x_{i}+\frac{h_{i+1}}{2},\, y_{i}+\frac{h_{i+1}}{2}f(x_{i}, \widehat{y}_{i})\right)= \widehat{y}_{i}+ h_{i+1}f_{i}+ \frac{h_{i+1}^2}{2}(f_x+ f_{i}\cdot f_y)+ O(h_{i+1}^3).[/math] Сравнивая последнее соотношение с (6.44), можно сделать вывод о том, что схема модифицированного метода Эйлера согласуется с разложением по формуле Тейлора вплоть до членов степени [math]h_{i+1}^2[/math] и поэтому является методом Рунге-Кутты второго порядка. 3. Коэффициенты семейства методов Рунге-Кутты (6.39),(6.40) удобно записывать в виде табл. 6.7. Табл. 6.8 соответствует классическому методу Рунге-Кутты четвертого порядка (6.41), а в табл. 6.9 приведены коэффициенты другого метода четвертого порядка (правило [math]3\!\!\not{\phantom{|}}\,8[/math]). Табл. 6.10 соответствует методу Рунге–Кутты (6.43) второго порядка, табл. 6.11 и 6.12 — методам третьего порядка. В табл. 6.13 даны коэффициенты метода Фельберга четвертого порядка, где формула для [math]\widetilde{y}_{i+1}[/math] имеет пятый порядок, а в табл. 6.14 — семистадийного метода Батчера шестого порядка. В табл. 6.15 приведены коэффициенты метода Рунге-Кутты-Мерсона четвертого порядка, где формула для [math]\widetilde{y}_{i+1}[/math] имеет пятый порядок для систем линейных дифференциальных уравнений с постоянными коэффициентами и третий порядок для нелинейных систем. Параметры метода Рунге-Кутты-Дормана-Принса пятого порядка приведены в табл. 6.16. Метод Фвльберга седьмого порядка и метод Дормана-Принса восьмого порядка изложены. 4. Для s-стадийных методов Рунге-Кутты справедливо следующее условие устойчивости, определяемое для задачи (6.12): [math]\left|1+ h\mu+ \frac{(h\mu)^2}{2}+ \frac{(h\mu)^3}{6}+ \frac{(h\mu)^4}{24}+ \ldots+\frac{(h\mu)^s}{s!}\right|<1.[/math] Если [math]\mu[/math] — действительное число, можно указать интервалы устойчивости: [math]h\mu\in (-2;0)[/math] при [math]s=1,\!2;~ h\mu\in(-2,\!51;0)[/math] при [math]s=3;~ h\mu\in (-2,\!78;0)[/math] при [math]s=4[/math] и т.д. Пример 6.6. Для задачи Коши [math]y'=x+y,~ y(0)=1,~ x\in[0;0,\!2][/math] выполнить два шага методом Рунге-Кутты второго порядка (6.43) и один (первый шаг) методом четвертого порядка (6.41). Шаг интефирования принять постоянным и равным 0,1. Сравнить оба численных решения с точным [math]y(x)= 2e^x-x-1[/math]. РешениеВ данной задаче [math]x_0=0,~ y_0=1[/math]. Применим схему (6.43). Первый шаг [math](i=0)[/math]. В соответствии с постановкой задачи принимается [math]h_1=0,\!1[/math]. Вычислим [math]\begin{aligned}K_{1,0}&= f(x_0,\widehat{y}_0)= f(x_0,y_0)= x_0+ y_0=0+1=1;\\[2pt] K_{2,0}&= f\! \left(x_0+ \frac{h_1}{2},\, y_0+\frac{h_1}{2}K_{1,0}\right)= 0,\!05+ 1,\!05= 1,\!1. \end{aligned}[/math] Тогда [math]\widehat{y}_1= \widehat{y}_0+ h_1K_{2,0}= 1+0,!1\cdot 1,!1= 1,\!11[/math]. Второй шаг [math](i=1)[/math]. Так как шаг постоянный, то [math]h_2=0,!1[/math], a [math]x_1=0,\!1[/math] и, следовательно, [math]\begin{aligned}K_{1,1}&= f(x_1,\widehat{y}_1)= x_1+ \widehat{y}_1= 0,\!1+1,\!11= 1,\!12\,;\\[2pt] K_{2,1}&= f\! \left(x_1+ \frac{h_2}{2},\, \widehat{y}_1+ \frac{h_2}{2}K_{1,1}\right)= 0,\!1+ \frac{0,\!1}{2}+ 1,\!11+ 0,\!05\cdot 1,\!12= 1,\!316\,. \end{aligned}[/math] Тогда [math]\widehat{y}_2= \widehat{y}_1+ h_2K_{2,1}= 1,\!11+ 0,\!1\cdot 1,\!316= 1,\!2416[/math]. Точное решение: [math]y(x_2)= 1,\!2428055[/math]. Отличие приближенного решения от точного составляет [math]0,\!096\%[/math]. Выполним теперь первый шаг методом четвертого порядка по формуле (6.41), записанной при [math]i=0\colon[/math] [math]\widehat{y}_1= \widehat{y}_0+ \frac{h_1}{6} (K_{1,0}+ 2K_{2,0}+ 2K_{3,0}+ K_{4,0}).[/math] При этом вычислим [math]K_{1,0}= x_0+ y_0=0+1=1[/math]; [math]\begin{aligned}K_{2,0}&= \left(x_0+ \frac{h_1}{1}\right)+ \left(\widehat{y}_0+ \frac{h_1}{2}\cdot K_{1,0}\right)= 0,\!05+ (1+0,\!05\cdot1)=1,\!1;\\[2pt] K_{3,0}&= \left(x_0+ \frac{h_1}{1}\right)+ \left(\widehat{y}_0+ \frac{h_1}{2}\cdot K_{2,0}\right)= 0,\!05+ (1+0,\!05\cdot1)=1,\!105;\\[2pt] K_{4,0}&= (x_0+ h_1)+ (\widehat{y}_0+ h_1\cdot K_{3,0})= 0,\!1+ (1+0,\!1\cdot1,\!105)=1,\!2105;\\[2pt] & \widehat{y}_1= 1+ \frac{0,\!1}{6} (1+ 2\cdot 1,\!1+ 2\cdot 1,\!105+ 1,\!2105)= 1,\!1103416. \end{aligned}[/math] Точное решение: [math]y(0,\!1)= 1,\!1103418[/math]. Сопоставляя точное решения [math]y(x_1)[/math] в точке [math]x_1[/math] с результатами [math]\widehat{y}_1[/math], полученными по схемам второго и четвертого порядка, можно увидеть, что для схемы четвертого порядка результат [math]\widehat{y}_1[/math] совпадает с точным решением в семи цифрах, а для схемы второго порядка — в трех знаках. В заключение данного раздела отметим преимущества и недостатки одношаговых и многошаговых схем: 1. Одношаговые схемы позволяют достаточно легко строить алгоритмы с переменным шагом и заданной точностью, что необходимо при расчете функций в областях их больших изменений. Однако они имеют недостаток — приближенные значения искомой функции, рассчитанные на предьщущих шагах, в этих схемах не используются. 2. В многошаговых схемах, напротив, учитываются значения искомой функции, рассчитанные на предыдущих шагах интегрирования, что повышает порядок точности этих схем без использования дополнительных точек на отрезке [math][x_{i},x_{i+1}][/math]. Однако эти схемы требуют усложнения расчетных алгоритмов при реализации переменного шага. Замечание. В явных одношаговых методах обеспечить заданную точность на каждом очередном шаге из текущей точки [math](x_{i}, \widehat{y}_{i})[/math] можно по следующей методике. 1. Задать начальную величину шага [math]h^{(0)}= x_{i+1}^{(0)}-x_{i}[/math] и рассчитать значение [math]\widehat{y}_{i+1}^{(0)}[/math] в точке [math]x_{i+1}^{(0)}[/math]. 2. Вычислить величину [math]\widehat{y}_{i+1}^{(1)}[/math] в точке [math]x_{i+1}^{(0)}[/math] с половинным значением шага [math]h^{(1)}= \frac{1}{2}h^{(0)}[/math] (до достижения точки [math]x_{i+1}^{(0)}[/math] выполняются два шага). 3. Если [math]\bigl|\widehat{y}_{i+1}^{(1)}-\widehat{y}_{i+1}^{(0)}\bigr|< \varepsilon[/math], то перейти к нахождению приближенного решения в точке [math]x_{i+2}^{(0)}= x_{i+1}^{(0)}+ h^{(0)}[/math]. Иначе осуществить еще одно дробление шага, т.е. вычислить [math]h^{(2)}= \frac{1}{2}h^{(1)}[/math], и найти [math]\widehat{y}_{i+1}^{(2)}[/math] в точке [math]x_{i+1}^{(1)}= x_{i}+ h^{(1)}[/math]. 4. Сравнить решения, полученные с шагом [math]h^{(1)}[/math] и [math]h^{(2)}[/math] в точке [math]x_{i+1}^{(1)}[/math]. Если желаемая точность достигнута, расчеты продолжить с шагом [math]h^{(1)}[/math], а если нет, то процесс дробления шага продолжить до достижения заданной точности. Принцип аналогийДанный метод следует из обобщения интегрально-интерполяционного принципа. Он основан на применении интерполяционных интегрально-дифференциальных сплайн-функций, преобразования которых — параметрические соотношения, аппроксимационные формулы для производных и квадратурные формулы — используются для получения разностных схем решения задачи Коши. Приведем некоторые из этих соотношений, соответствующих задаче аппроксимации некоторой произвольной функции [math]\Phi(x)[/math] многочленами второй степени: а) параметрическое соотношение, связывающее интерполируемую функцию [math]\Phi(x_{i})[/math] и ее производную [math]\overline{m}_{i}= \Pho'(x_{i})\colon[/math] [math]\frac{\overline{m}_{i}}{h_{i+1}}-\frac{\overline{m}_{i-1}}{h_{i}}= \frac{\Delta \Phi_{i}}{h_{i+1}^2}-\frac{\Delta \Phi_{i-1}}{h_{i}^2}\,;[/math] (6.46) б) формула аппроксимации первой производной в точке [math]x_{i}[/math] на трехточечном нерегулярном шаблоне [math](x_{i-1}, x_{i},x_{i+1})\colon[/math] [math]\overline{m}_{i+1}= \frac{}{}\! \left[\delta_{i+1}\Phi_{i-1}- \frac{(H_{i}^{i})^2}{\Pi_{i}^{i+1}}\Phi_{i}+ \frac{H_{i}^{2(i+1)}}{h_{i+1}}\Phi_{i-1}\right]~ \left(\frac{h_{i}^2}{6} \delta_{i+1} (1+\delta_{i+1})M_{3,i}\right)\!.[/math] (6.47) где [math]H_{si}^{p(i+1)}= s\cdot h_{i}+ p\cdot h_{i+1};~ s,\,p[/math] — целые числа [math]\left(H_{si}^{p(i+1)}= h_{i}(s+ p\cdot \delta_{i+1}),~ \delta_{i+1}= \frac{h_{i+1}}{h_{i}}\right)[/math]; в) квадратурная формула экстраполяционного типа, имеющая четвертый порядок: [math]I_{i}^{i+1}= \frac{h_{i+1}^2}{6H_{i-1}^{i}} \! \left[\frac{H_{3i}^{2(i+1)}}{h_{i-1}}\Phi_{i-2}-\frac{3(H_{i-1}^{i})^2+2h_{i+1}H_{i-1}^{i}}{\Pi_{i-1}^{i}}\Phi_{i-1}+ \frac{6h_{i}H_{i-1}^{i}+ h_{i+1}(3h_{i-1}+ 6h_{i}+ 2h_{i+1})}{\Pi_{i}^{i+1}}\Phi_{i}\right]\!.[/math] (6.48) Эта формула при [math]h_{i+1}= h=\text{const}[/math] упрощается: [math]I_{i}^{i+1}= \frac{h}{12} (5\Phi_{i-2}-16\Phi_{i-1}+ 23\Phi_{i})\quad \left(\frac{3}{8} h^4M_{3,i}\right)\!;[/math] г) двухинтервальная (трехточечная) обобщенная на нерегулярный шаблон квадратурная формула парабол четвертого порядка на отрезке [math][x_{i-1},x_{i+1}][/math], выражающая сумму "удельных" значений интегралов через линейную комбинацию значений [math]\Phi_{i-1}, \Phi_{i}, \Phi_{i+1}\colon[/math] [math]3\! \left(\frac{1}{h_{i}^2}I_{i-1}^{i}+ \frac{1}{h_{i+1}^2} I_{i}^{i+1}\right)= \frac{1}{h_{i}}\Phi_{i-1}+ 2\! \left(\frac{1}{h_{i}}+ \frac{1}{h_{i+1}}\right)\!\Phi_{i}+ \frac{\Phi_{i+1}}{h_{i+1}}\,;[/math] (6.49) д) одноинтервальная (трехточечная) функциональная квадратурная формула четвертого порядка на отрезке [math][x_{i-1},x_{i+1}][/math] [math]I_{i}^{i+1}= \frac{h_{i+1}^3}{6H_{i}^{i+1}}\! \left[-\frac{1}{h_{i}}\Phi_{i-1}+ \frac{H_{i}^{i+1}H_{3i}^{i+1}}{\Pi_{i}^{i+1} h_{i+1}}\Phi_{i}+ \frac{H_{3i}^{2(i+ 1)}}{h_{i+1}^2}\Phi_{i+1}\right]~ \left(\frac{h_{i}^2h_{i+1}^2}{72} \delta_{n+1} (\delta_{n+1}+2)M_{3,i}\right)\!.[/math] (6.50) которая при [math]h_{i+1}= h=\text{const}[/math] преобразуется к более простому виду: [math]I_{i}^{i+1}= \frac{h}{12} \bigl(-\Phi_{i-1}+ 8\Phi_{i}+ 5\Phi_{i+1}\bigr).[/math] Анализ остаточного слагаемого, получающегося при оценке погрешности квадратурной формулы (6.50) (см. мажоранту в скобках справа от данной формулы), указывает на то, что при сгущающейся вправо сетке (условная аппроксимация интеграла) порядок точности формулы может быть повышен на единицу без увеличения количества точек шаблона. Действительно, при [math]\delta_{i+1}< h_{i}[/math] легко получить приближенное условие повышения порядка: [math]h_{i+1} \leqslant h_{i}^2[/math]. Таким образом, в данном случае порядок повышается за счет значительного уменьшения шага, такого, что [math]h_{i+1}[/math] не превышает квадрата шага [math]h_{i}[/math]. Этот факт может быть использован при решении задачи Коши на нерегулярном шаблоне; е) две одноинтервальные (двухточечные) квадратурные формулы функционально-дифференциального типа: [math]I_{i}^{i+1}= \frac{h_{i+1}}{3} (\Phi_{i}+ 2\Phi_{i+1})-\frac{h_{i+1}^2}{6}\Phi'_{i+1}\quad \left(\frac{5}{24} h_{i}^4M_{3,i}\right)\!,[/math] (6.51) [math]I_{i}^{i+1}= h_{i+1}\Phi_{i}+ \frac{h_{i+1}^2}{6} (2\Phi'_{i}+ \Phi'_{i+1})\quad \left(\frac{1}{24}h_{i}^4M_{3,i}\right)\!;[/math] (6.52) ж) интегрально-функциональное параметрическое соотношение, осуществляющее связь интегралов и производных: [math]\frac{I_{i}^{i+1}}{h_{i+1}}-\frac{I_{i-1}^{i}}{h_{i}}= \frac{1}{6} \bigl(h_{i+1} \overline{m}_{i+1}+ 2H_{i}^{i+1}\cdot \overline{m}_{i}+ h_{i}\cdot \overline{m}_{i-1}\bigr).[/math] (6.53) Дискретные явные или неявные схемы для решения задачи (6.9) получаются из приведенных соотношений либо посредством аппроксимации производной, либо путем соответствующей замены одних параметров, относящихся к параметрическому соотношению или к квадратурной формуле, на другие параметры. Этот способ построения схем называется способом аналогий. Так, заменяя в (6.46) [math]\Delta\Phi_{k}[/math] на [math]\Delta y_{k}[/math], а [math]\overline{m}_{k}[/math] на [math]f_{k}~ (k=i,i-1)[/math] и разрешая относительно [math]\widehat{y}_{i+1}[/math], получаем двухшаговую явную схему второго порядка (2Я2А): [math]\widehat{y}_{i+1}= \widehat{y}_{i}+ \delta_{i+1}^2 \Delta \widehat{y}_{i-1}+ h_{i+1}^2\! \left(\frac{f_{i}}{h_{i+1}}-\frac{f_{i-1}}{h_{i}}\right)\!,[/math] (6.54) которая при [math]h_{i+1}= h= \text{const}[/math] преобразуется к виду [math]\widehat{y}_{i+1}= 2\widehat{y}_{i}-\widehat{y}_{i-1}+ h(f_{i}-f_{i-1})[/math]. Комбинация схем (6.54) и (6.24) дает обобщенную на нерегулярный шаблон схему Адамса—Бэшфорта (2Я2В): [math]\widehat{y}_{i+1}= \widehat{y}_{i}+ \frac{h_{i+1}^2}{2}\! \left(\frac{H_{2i}^{i+1}}{\Pi_{i}^{i+1}} f_{i}-\frac{1}{h_{i}}f_{i-1}\right)\!,[/math] (6.55) которая при [math]h_{i+1}= h= \text{const}[/math] преобразуется к виду (6.32): [math]\widehat{y}_{i+1}= \widehat{y}_{i}+ \frac{h}{2}(3f_{i}-f_{i-1})[/math]. Аппроксимационная формула (6.47) определяет неявную двухшаговую схему второго порядка (2НЯ2): [math]\widehat{y}_{i+1}= \frac{h_{i+1}}{H_{i}^{2(i+1)}}\! \left[\frac{(H_{i}^{i+ 1})^2}{\Pi_{i}^{i+1}} \widehat{y}_{i}-\delta_{i+1} \widehat{y}_{i-1}+ H_{i}^{i+1} f(x_{i+1}, \widehat{y}_{i+1})\right]\!,[/math] (6.56) которая при [math]h_{i+1}= h= \text{const}[/math] принимает вид [math]\widehat{y}_{i+1}=-\frac{1}{3} \widehat{y}_{i-1}+ \frac{4}{3} \widehat{y}_{i}+ \frac{2}{3} h\, f(x_{i+1},\widehat{y}_{i+1}).[/math] (6.57) Замечание. Схема (6.57) является одним из методов Гира, получаемых на основе формул дифференцирования назад, связывающих значение производной в точке [math]x_{i+1}[/math] со значениями функции в предыдущих точках [math]x_{i},\ldots, x_{i-k+1}[/math] (см. рис.6.2). Схема (6.57) является схемой второго порядка, а неявный метод Эйлера — схемой первого порядка. Приведем более точные многошаговые схемы, записанные на регулярных шаблонах: – схема третьего порядка: [math]\widehat{y}_{i+1}= \frac{18}{11} \widehat{y}_{i}-\frac{9}{11} \widehat{y}_{i-1}+ \frac{2}{11} \widehat{y}_{i-2}+ \frac{6}{11} h\, f(x_{i+1}, \widehat{y}_{i+1});[/math] (6.58) – схема четвертого порядка: [math]\widehat{y}_{i+1}= \frac{48}{25} \widehat{y}_{i}-\frac{16}{25} \widehat{y}_{i-1}+ \frac{16}{25} \widehat{y}_{i-2}-\frac{3}{25}\widehat{y}_{i-3}+ \frac{12}{25}h\, f(x_{i+1}, \widehat{y}_{i+1});[/math] – схема пятого порядка: [math]\widehat{y}_{i+1}= \frac{300}{137} \widehat{y}_{i}-\frac{300}{137} \widehat{y}_{i-1}+ \frac{200}{137} \widehat{y}_{i-2}-\frac{75}{137}\widehat{y}_{i-3}+ \frac{12}{137} \widehat{y}_{i-4}+ \frac{60}{137}h\, f(x_{i+1}, \widehat{y}_{i+1});[/math] – схема шестого порядка: [math]\widehat{y}_{i+1}= \frac{360}{147} \widehat{y}_{i}-\frac{450}{147} \widehat{y}_{i-1}+ \frac{400}{147} \widehat{y}_{i-2}-\frac{225}{147}\widehat{y}_{i-3}+ \frac{17}{147} \widehat{y}_{i-4}-\frac{10}{147} \widehat{y}_{i-5}+ \frac{60}{147}h\, f(x_{i+1}, \widehat{y}_{i+1}).[/math] Известно, что схема (6.58) является A(α)-устойчивой с [math]\alpha\cong 68^{\circ}[/math]. Все рассматриваемые ниже схемы получаются из квадратурных формул (6.48)-(6.52), в которых вместо разности первообразных [math]F_{i+1}-F_{i}= I_{i}^{i+1}[/math] подставляется [math]\Delta \widehat{y}_{i}[/math], а вместо [math]\Phi_{k}(k=i-2,i-1,i)[/math] и [math]\Phi'_{k}[/math] подставляются [math]f_{k}[/math] и [math]f'_{k}[/math]. При этом порядки точности схем на отрезке [math][a,b][/math] соответствуют порядкам точности вычисления интегралов на всем отрезке интегрирования. Преобразуя таким образом квадратурную формулу (6.48), получаем явную трехшаговую схему третьего порядка (ЗЯЗ): [math]\widehat{y}_{i+1}= \widehat{y}_{i}+ \frac{h_{i+1}^2}{6H_{i-1}^{i}}\! \left[\frac{H_{3i}^{i+1}}{h_{n-1}} f_{i-2}-\frac{3(H_{i-1}^{i})^2+ 2h_{i+1} H_{i-1}^{i}}{\Pi_{i-1}^{i}}f_{i-1}+ \frac{6h_{i}H_{i-1}^{i}+ h_{i+1} (H_{3(i-1)}^{6i}+ 2h_{i+1})}{\Pi_{i}^{i+1}}f_{i}\right]\!.[/math] (6.59) При [math]h_{i+1}= h= \text{const}[/math] она является традиционной, получающейся интегрально-интерполяционным методом (6.33): [math]\widehat{y}_{i+1}= \widehat{y}_{i}+ \frac{h}{12}\bigl(5f_{i-2}-16 f_{i-1}+ 23f_{i}\bigr).[/math] Из обобщенной формулы (6.49) следует первая двухшаговая неявная схема третьего порядка (2НЯЗА), которая при [math]h_{i+1}= h= \text{const}[/math] преобразуется к виду (6.29) [math]\widehat{y}_{i+1}= \widehat{y}_{i}-\delta_{i+1}^2 \Delta \widehat{y}_{i-1}+ \frac{h_{i+1}^2}{3}\! \left[\frac{f_{i-1}}{h_{i}}+ 2\! \left(\frac{1}{h_{i}}+ \frac{1}{h_{i+1}}\right)\! f_{i}+ \frac{1}{h_{i+1}} f(x_{i+1}, \widehat{y}_{i+1})\right]\!.[/math] (6.60) Одноинтервальная (трехточечная) квадратурная формула (6.50) определяет вторую двухшаговую неявную схему третьего порядка (2НЯЗБ) [math]\widehat{y}_{i+1}= \widehat{y}_{i}+ \frac{h_{i+1}^3}{6H_{i}^{i+1}}\! \left[-\frac{f_{i-1}}{h_{i}}+ \frac{H_{i}^{i+1} H_{3i}^{i+1}}{\Pi_{i}^{i+1} h_{i+1}}f_{i}+ \frac{H_{3i}^{2(i+1)}}{h_{i+1}^2}f(x_{i+1}, \widehat{y}_{i+1})\right]\!.[/math] (6.61) Эта схема в соответствии с вышеприведенной оценкой точности интеграла при [math]\delta_{i+1} \leqslant h_{i}[/math] имеет не третий, а четвертый порядок точности. При [math]h=\text{const}[/math] эта схема преобразуется к виду (6.36): [math]\widehat{y}_{i+1}= \widehat{y}_{i}+ \frac{h}{12} \bigl(-f_{i-1}+ 8f_{i}+ 5f(x_{i+1}, \widehat{y}_{i+1})\bigr).[/math] И наконец, последние приводимые здесь схемы следуют из квадратурных формул (6.51), (6.52) (неявные одношаговые схемы) и параметрического соотношения (6.53) (неявная двухшаговая схема): [math]\begin{aligned}\widehat{y}_{i+1}&= \widehat{y}_{i}+ \frac{h_{i+1}}{3} \bigl[f_{i}+ 2f(x_{i+1}, \widehat{y}_{i+1})\bigr]-\frac{h_{i+1}^2}{6}f'(x_{i+1}, \widehat{y}_{i+1});\\[2pt] \widehat{y}_{i+1}&= \widehat{y}_{i}+ h_{i+1}f_{i}+ \frac{h_{i+1}^2}{6} \bigl[2f'_{i}+ f'(x_{i+1}, \widehat{y}_{i+1})\bigr];\\[2pt] \widehat{y}_{i+1}&= \widehat{y}_{i}+ \delta_{i+1} \Delta \widehat{y}_{i+1}+ \frac{h_{i+1}}{6} \bigl[h_{i} f'_{i-1}+ 2H_{i}^{i+1}f'_{i}+ h_{i+1} f'(x_{i+1}, \widehat{y}_{i+1})\bigr], \end{aligned}[/math] где [math]f'=f_{x}+f_{y}\cdot f[/math]. Эти схемы обозначаются 1НЯЗА, 1НЯЗБ, 2НЯЗД соответственно. При [math]h_{i+1}= h= \text{const}[/math] последняя схема преобразуется к виду [math]\widehat{y}_{i+1}=-\widehat{y}_{i-1}+ 2\widehat{y}_{i}+ \frac{h^2}{6} \bigl[f'_{i-1}+ 4f'_{i}+ f'(x_{i+1}, \widehat{y}_{i+1})\bigr].[/math] Для применения разностных схем, полученных с помощью четырех различных подходов, можно использовать общую методику. Методика решения задачи Коши методами сеток1. В соответствии с требуемой точностью решения задачи Коши выбрать численный метод (схему) определенного порядка аппроксимации и типа (явный или неявный, одношаговый или многошаговый). 2. Задать сетку [math]\Omega_{n}[/math] или выбрать ее параметр ([math]h[/math], если сетка равномерная; [math]h_{i+1},~ i=\overline{0,n-1}[/math], если сетка неравномерная). Если выбранный метод является ограниченно устойчивым, то для задания параметра сетки оценить величину критического шага. 3. Если выбранный метод многошаговый, найти требуемое число "разгонных" точек. 4. Найти приближенное решение задачи Коши во всех точках сетки [math]\Omega_{n}[/math]. Если выбран явный метод, каждая следующая точка [math]\widehat{y}_{i+1}[/math] рассчитывается по явной формуле, а если неявный, то с помощью решения в общем случае нелинейного алгебраического уравнения (6.11) одним из численных методов. 5. Оценить достигнутую точность. В случае невыполнения требований повторить расчеты, меняя значение шага или используя более точный метод.
Математический форум (помощь с решением задач, обсуждение вопросов по математике).
Если заметили ошибку, опечатку или есть предложения, напишите в комментариях.
|
Часовой пояс: UTC + 3 часа [ Летнее время ] |