Численное дифференцирование и интегрирование

Тема 4. Численное ДИФФЕРЕНЦИРОВАНИЕ И ИНТЕГРИРОВАНИЕ 1. Основные формулы численного дифференцирования 1.1. Трехточечная формула для у’(x) при произвольном х , (4.1) где . 1.2. Трехточечные формулы для у’(x) в узлах интерполяции ; ; (4.2) . 1.3. Трехточечная формула для у’’(x) , где . 1.4. Частные производные ; ; ; ; = =. Для всех приведенных формул частных производных погрешность составляет О(h2) (O(ph) для смешанной производной; р – шаг по переменной у). 1.5. Обусловленность численного дифференцирования. Поскольку при численном дифференцировании вычитаются два близких числа, то данная задача плохо обусловлена. Пусть | y”’()|  М3 и абсолютная погрешность вычисления у(х) не превосходит Е. Тогда при использовании формулы (4.2) для вычисления у’(x1) справедлива оценка . Величина обусловленности представляет собой коэффициент при Е, следовательно,  = 1/ h. Для улучшения обусловленности можно увеличить h, используя формулу оптимального выбора шага h. При использовании формулы (4.2) для вычисления у’(x1) справедливы оценки: ; ЧИСЛ. ДИФФ.  1.04. (4.3) Пример 4.1. Вычислим значение производной функции f(x) = ex по формуле (4.2) в точке х = 1 при различных значениях h на ЭВМ, вычисляющей с точностью до 6 значащих цифр. Результаты – в таблице 4.1. Таблица 4.1 Численное дифференцирование у = ех при различных h h 0.1 0.05 0.025 0.01 0.005 0.001 0.0001 y’(x) 2.72285 2.7194 2.7186 2.7185 2.719 2.72 2.7 (y’*) 4.5710–3 1.1210–3 3.1810–4 2.1810–4 7.1810–4 1.7210–3 –0.018 На рис. 4.1 данные таблицы приведены в графическом виде с использованием логарифмического масштаба по обеим осям (Ох – lg(h), Oy – lg()). Из рисунка видно, что минимальная погрешность имеет место при h  0.01 (lg h = –2). Вычисления по формулам (4.3) дают близкий результат: = 0.01033. Пример 4.2. Вычислить значение первой производной функции у(х) = ln x по формуле (4.2) в точке х = 1.2 с оптимальным выбором шага. Положить Е = 10 –11. Для вычисления оценки y”’() воспользуемся формулой . Построим таблицу разделенных разностей у(х) при хi = 1; 1.2; 1.4; 1.6. Получим 3 = 0.1573. Отсюда h opt = = 3.167722  10 –4 ; ЧИСЛ. ДИФФ.  4.335  10 –8. По второй из формул (4.2) получаем y’(1.2) . Сравнивая с реальным значением , приходим к выводу, что погрешность численного дифференцирования не превзошла оценочную. 2. Простые формулы численного интегрирования 2.1. Постановка задачи численного интегрирования. В прикладных исследованиях часто возникает необходимость вычисления значения определенного интеграла , который может выражать площадь, объем, работу переменной силы и т.д. Если функция f(x) непрерывна на [a, b] и ее первообразную F(x) удается выразить через известные функции, то для вычисления интеграла можно воспользоваться формулой Ньютона-Лейбница: I = F(b) – F(a). (4.4) К сожалению, в подавляющем большинстве случаев получить значение интеграла с помощью формулы (4.4) или других аналитических методов не удается. Так, интеграл широко используется при исследовании процессов массо- и теплообмена, в статистической физике и теории вероятностей. Однако его значение не может быть выражено в виде конечной комбинации элементарных функций. При разработке численных методов нахождения значения определенных интегралов наиболее распространенным подходом является аппроксимация. Подынтегральную функцию f(x) на отрезке [a, b] приближают некоторой функцией Y(x), которая легко интегрируется аналитически. Затем полагают . Если функция f(x) задана аналитически, то ставится вопрос об оценке погрешности такой замены. 2.2. Формула трапеции. Пусть функция f(x) интерполируется на [a, b] многочленом 1-й степени: Y(x) = P1 (x) = y0 + (х – x0) 1[x0, x1], где x0 = a, x1= b. Обозначим h = x1– x0. После интегрирования и преобразования получим следующую формулу:  [f(a) + f(b)] = [f(x0) + f(x1)] = I*, (4.5) называемой формулой трапеций. Геометрически она означает замену площади криволинейной трапеции обычной трапецией (рис. 4.2), отсюда название формулы. Воспользовавшись известной оценкой погрешности интерполирования многочленом P1(x) несложно получить оценку погрешности формулы (4.5): RТРАП = I – I* = , [ x0 , x1]. (4.6) Знак “–” означает, что при f’’ > 0 формула (4.5) дает значение интеграла с избытком, а при f ’’ < 0 – с недостатком. На рис. 4.1 функция f(x) вогнута, т.е. f ’’ < 0, следовательно. I* < I. При f ’’ = 0, т.е. если f(x) = P1(x), формула (4.5) не содержит погрешности и абсолютно точна. 2.3.Формула Симпсона. Пусть функция f(x) интерполируется на [a, b] многочленом 2-й степени: Y(x) = P1 (x) = y0 + (х – x0) ( 1[x0, x1] + (х – x1) 2[x0, x1, x2] ), где x0 = a, x1 = (a + b) / 2, x2 = b, h = xi+1 – xi. После интегрирования и преобразования получим следующую формулу:  [ f (x0) + 4f (x1) + f (x2) ] = I*, (4.7) называемой формулой Симпсона. Ее геометрический смысл представлен на рис. 4.3. Оценка погрешности формулы (4.7) имеет вид: RСИМП = I – I* = , [ x0 , x2]. (4.8) Знак “–” означает, что при f (4) > 0 формула (4.7) дает значение интеграла с избытком, а при f (4) < 0 – с недостатком. На рис. 4.4 функция f(x) представляет собой многочлен 4-й степени с постоянным значением f (4) < 0. Из рисунка видно, что площадь фигуры под графиком f(x) (окрашенная область) больше, чем под Y(x). Следовательно, I* < I. При f (4) = 0, т.е. если f(x) = P3(x), формула (4.6) не содержит погрешности и абсолютно точна. Это видно на рис. 4.3 – площади под f(x) и под Y(x) абсолютно одинаков, хотя сами функции не совпадают. 3. Составные формулы численного интегрирования Описанные формулы численного интегрирования легко обобщить на случай применения многочленов более высоких степеней и тем самым добиться улучшения точности решения задачи (подобный метод широко применяют, при этом, в частности, получаются так называемые квадратурные формулы Ньютона-Котеса). Однако, более перспективен другой способ, основанный на кусочной аппроксимации подынтегральной функции многочленами невысокой степени. 3.1. Кусочно-гладкая аппроксимация функций. Пусть функция интерполируется на отрезке [a, b]. Метод решения этой задачи с помощью единого для всего отрезка многочлена Рn(х) называют глобальной многочленной интерполяцией. При всей привлекательности перспективы иметь такой единый многочлен, существуют весьма веские причины, по которым глобальная интерполяция многочленами высокой степени в вычислительной практике, как правило, не используется. Одна из основных причин этого заключается в проблеме сходимости аппроксимации при росте n. Рассмотрим пример. Пример 4.3. Непрерывную функцию у = | x | надо приблизить на отрезке [–1, 1] с высокой точностью. На рис. 4.5 и 4.6. приведены графики интерполяционных многочленов 5 и 13 степеней, соответственно. Видно, что при сравнительно невысокой степени n = 5 велика погрешность аппроксимации в точке х = 0. При более высоких степенях, начиная с n = 9, появляются “паразитические” всплески вблизи точек х = 1, нарастающие при увеличении n. В результате высокая точность, порядка  = 0.001, оказывается недостижимой. В известной степени точность глобальной многочленной интерполяции можно повысить за счет оптимального расположения узлов, однако полностью этим проблема не снимется. Вторая проблема заключается в ухудшении обусловленности задачи интерполяции при росте n, и соответственно, в повышении чувствительности к погрешности округления. В значительной мере проблемы глобальной аппроксимации могут быть сняты применением кусочно-гладкой, в частности – кусочно-полиномиальной интерполяции. Суть подхода заключается в следующем. Исходный отрезок [a, b] разбивается на несколько отрезков меньшей длины, на каждом из которых функция приближается своей гладкой кривой, например, многочленом. Низкая степень многочленов гарантирует отсутствие нежелательных эффектов, подобных рассмотренным в примере 4.3, а малая величина каждого отрезка способствует увеличению точности аппроксимации. Так, в примере 4.3 исходную функцию можно абсолютно точно приблизить двумя многочленами 1-й степени на отрезках [–1, 0] и [0, 1], соответственно. Замечание. Подобный подход также имеет недостатки, главным из которых является негладкость интерполяционных кривых в точках “стыковки” отдельных многочленов. Данный недостаток может быть преодолен использованием так называемой сплайн-аппроксимации, при которой параметры многочленов подбираются из условия гладкости сочленения (совпадения первых производных) отдельных звеньев (сплайнов) в точках “склеивания”. Подобная интерполяция используется при решении многих задач, в частности, в компьютерной графике и машиностроении при построении линий сложной формы. На рис. 4.7 приведены графики различных видов интерполяции. Для нужд численного интегрирования негладкость кусочной аппроксимации, вообще говоря, помехой не является, поэтому методы аппроксимации сплайнами подробно нами рассматриваться не будут. 3.2. Составная формула трапеции. Следуя идее кусочно-гладкой аппроксимации, разобьем отрезок [a, b] на n равных частей точками x0 = a, xn = b, h = xi+1 – xi, обозначим yi = f(xi), и применим формулу (4.5) к каждому отрезку [xi, xi+1]. В результате получим составную формулу трапеций: I* = h . (4.9) Ее погрешность выражается формулой = I – I* = , [a , b]. (4.10) 3.3. Составная формула Симпсона. Разобьем отрезок [a, b] на четное число равных частей точками x0 = a, xn = b, h = xi+1 – xi, обозначим yi = f(xi), и применим формулу (4.7) к каждому отрезку [xi, xi+2]. В результате получим составную формулу Симпсона: I* = [ (y0 + yn) + 4 (y1 + y3 + … yn–1) + 2 (y2 + y4 + … yn–2) ]. (4.11) Ее погрешность выражается формулой = I – I* = , [a , b]. (4.12) 4. Метод дробления шага. Правило Рунге Формулы (4.10) и (4.12) для практики неудобны, т.к. значения у(m) () обычно неизвестны. Правило Рунге позволяет найти достаточно точные оценки погрешности, используя значения I*, вычисленные при различных h. Пусть m – порядок точности применяемого метода численного интегрирования. Напомним, что для формулы трапеций m = 2, а для формулы Симпсона m = 4. Пусть далее, Ih и Ih/2 – значения интеграла, найденные при величине шага h и h/2, соответственно. Имеет место формула практической оценки погрешности согласно правилу Рунге: . (4.13) Если значение правой части приближенного равенства окажется меньше предельно допустимой погрешности , то можно считать Ih/2 ответом задачи. В противном случае – требуемая точность не достигнута и необходимо использовать более мелкий шаг, т.е. продолжить процесс дробления шага. Пример 4.4. Вычислить с точностью  = 10–5 значение интеграла , используя правило Рунге. 1) Формула трапеций. Для определения начального шага воспользуемся формулой (4.10). При заданной максимальной абсолютной погрешности  значение h можно найти по формуле . Величину у”() можно оценить по известной формуле , построив таблицу разделенных разностей, либо положить его равным какому-нибудь близкому по порядку числу, например 1. В любом случае мы угадаем лишь порядок оптимального значения h. Итак, полагая у”() = 1, получим h = 8.16510–3 . Отсюда число разбиений n = 220. По формуле (4.10) получим I h = 2.5651958. Уменьшаем шаг вдвое – I h/2 = 2.5650110. Отсюда = –6.2  10 –5. Значит, при таком разбиении требуемая точность не достигается, надо продолжить процесс дробления шага. Полностью процесс вычислений приведен в таблице 4.2. Окончательный результат I  2.5649532 получен при n = 1760. Сравнивая с истинным значением I = 2.56494936, убеждаемся, что применение правила Рунге позволило получить результат с требуемой точностью. Таблица 4.2 Правило Рунге для формулы трапеций n I* I – Ih/2 220 2.5651958 440 2.5650110 –6.2  10 –5 880 2.5649548 –1.6  10 –5 1760 2.5649532 –3.9  10 –6 2) Формула Симпсона. Начальное значение h равно Значит, отрезок [–0.8, 1] надо разбить на 10 частей. По формуле (4.11) получим I h = 2.586503286. Уменьшаем шаг вдвое – I h/2 = 2.56755997441. Отсюда = –1.26  10 –3. Значит, при таком разбиении требуемая точность не достигается, надо продолжить процесс дробления шага. Полностью процесс вычислений приведен в таблице 4.3. Окончательный результат I  2.5649504 получен при n = 160. Сравнивая с истинным значением I = 2.56494936, убеждаемся, что применение правила Рунге позволило получить результат с требуемой точностью. Таблица 4.3 Правило Рунге для формулы Симпсона n I* I – Ih/2 10 2.5865033 20 2.5675600 –1.26  10 –3 40 2.5651762 –1.6  10 –4 80 2.5649654 –1.4  10 –5 160 2.5649504 –1.0  10 –6 Замечание 1. В обоих случаях понадобилось повторить процесс дробления несколько раз. Такая медленная сходимость объясняется наличием особенности подынтегральной функции вблизи отрезка интегрирования в точке х = –0.95, в связи с чем, значения у(m) () растут вместе с уменьшением h. При интегрировании более гладких функций сходимость обычно бывает значительно лучше. Замечание 2. При приближенном вычислении интегралов, абсолютная величина | I | которых мала, использование значения абсолютной погрешности  для контроля точности необоснованно, более приемлема относительная погрешность . Поэтому при использовании правила Рунге рекомендуется следующий критерий останова процесса дробления шага: если | Ih/2 |  1, то заканчивать при , если | Ih/2 | < 1, то заканчивать при . 5. Тестовые примеры для программ численного интегрирования Приведем данные трассировки программ вычисления интеграла примера 4.4 при n = 8. При таком разбиении используются следующие узловые точки: Таблица 4.4 Узловые точки численного интегрирования примера 4.4 х f(x) x f(x) x f(x) –0.8 6.6666667 –0.125 1.2121212 0.55 0.6666667 –0.575 2.6666667 0.1 0.9523810 0.775 0.5797101 –0.35 1.6666667 0.325 0.7843137 1 0.5128205 Значение Ih, вычисленное при таком разбиении по формуле трапеции составит 2.7266106655, причем у0 + уn = 7.1794871795, = 8.5285260349. При дроблении шага получим Ih/2 = 2.6093879352. При использовании формулы Симпсона получим следующие результаты: Ih = 2.6041622061; (y1 + y3 + … yn–1) = 5.2428117492; (y2 + y4 + … yn–2) = 3.2857142857; Ih/2 = 2.5703136918. Тема 5. Численное решение ОБЫКНОВЕННЫХ дифференциальных уравнений (ОДУ) Решение дифференциальных уравнений – одна из наиболее часто встречающихся вычислительных задач, которая возникает как непосредственно при моделировании реальных объектов, так и в качестве промежуточной задачи. Как правило, точное решение не удается выразить через элементарные функции в конечном виде. Если провести грубое сравнение доли “берущихся” интегралов среди общего числа существующих типов интегралов и доли дифференциальных уравнений, решаемых аналитически среди всех известных типов дифференциальных уравнений, то окажется, что для уравнений эта доля намного меньше, чем для интегралов. Иными словами, задача численного решения дифференциальных уравнений весьма и весьма актуальна. 1. Классификация численных методов решения дифференциальных уравнений Дифференциальное уравнение и его численное решение. Напомним, что дифференциальным уравнением называется уравнение вида F(x, y, y’, y’’,…) = 0, (5.1) а его решением – такая функция у = у(х), которая при ее подстановке в F обращает уравнение (5.1) в тождество. Максимальный порядок производной в выражении F называется порядком уравнения. Уравнение (5.1) не имеет однозначного решения. Чтобы однозначно определить решение, требуются дополнительные условия. По типу таких условий различают одноточечные задачи (задачи с начальными условиями или задача Коши) и многоточечные. В первом случае все дополнительные условия задаются в одной, обычно начальной, точке. Пример 5.1. Дано уравнение 2-го порядка с двумя начальными условиями В многоточечных задачах дополнительные условия задаются в нескольких точках. Чаще всего – в граничных, тем самым ставится так называемая краевая задача. Пример 5.2. Дано уравнение 2-го поряка с условиями на концах отрезка определения решения – краевая задача Задачи обоих примеров имеют одинаковое аналитическое решение у = 1.125 e2x – 0.125 e–2x – 2x. Среди множества задач, связанных с решением дифференциальных уравнений наиболее важной является задача Коши для уравнения 1-го порядка, к ней сводятся многие другие задачи. Поэтому, именно ей будет уделено основное внимание в настоящей теме. Если аналитическое решение дифференциального уравнения представляет собой некоторое выражение, явно или неявно задающее искомую функцию у(х), то численное решение представляет собой таблицу, определяющую значение у(х) в некоторых заданных узлах хi. В дальнейшем построенная таблица может быть использована для получения решения в промежуточных точках, например, с применением аппарата аппроксимации функций. В таблице 5.1 приведено численное решение задачи примеров 5.1 и 5.2 в равноотстоящих точках. Таблица 5.1. Численное решение примера 5.1 xi yi y’i 1 0.5 0.25 1.2790 1.8613 0.50 2.0121 4.2081 0.75 3.5140 8.1396 1.00 6.2958 14.6592 1.25 11.1950 25.4311 1.50 19.5900 43.2049 1.75 33.7511 72.5173 2.00 57.4206 120.8503 Численные методы. Большинство численных методов решения задачи Коши для уравнения 1-го порядка, сводятся к процедуре вычисления уi+1 = y (x i + h) при известных значениях уi = y (x i ) и некоторых других параметров задачи. Будем считать, что уравнение (5.1) разрешено относительно у’ и мы имеем задачу Коши вида . (5.2) Зная начальную точку (х0 , у0 ), можно с помощью основной процедуры численного метода определить следующую точку (х1 , у1 ) решения. А дальше вычисления повторяются – положив (х1 , у1 ) новой начальной точкой мы можем тем же способом найти (х2 , у2 ) и так продолжить до конца отрезка задания решения. Основой для получения расчетных формул того или иного численного метода является эквивалентное представление уравнения (5.2) при начальной точке (хi , уi ) в интегральном виде: , (5.3) который получается интегрированием обеих частей уравнения (5.2). Для построения метода надо подобрать подходящую аппроксимацию интеграла в (5.3) с помощью какой-либо формулы численного интегрирования. Так как большинство квадратурных формул интерполяционного типа сводятся к вычислению значений у(х) в специально подобранныx узлах, то при аппроксимации интеграла вместо (5.3) будем иметь приближенное равенство уi+1 – уi  Ф(хi , h, уi+1–k ,… уi, уi+1), (5.4) где Ф – аппроксимация интеграла в (5.3). Теперь задача сводится к тому, чтобы выразить из алгебраического уравнения (5.4) искомую величину уi+1. Тип численного метода в основном определяется способом построения функции Ф и процедурой нахождения уi+1 из (5.4). При аппроксимации интеграла (5.3) в общем случае используются найденные ранее k значений искомой функции: уi+1–k,…, уi. Поэтому такие методы называются k-шаговыми. При k = 1 уравнение (5.4) упрощается и приобретает вид уi+1 – уi  Ф(хi , h, уi, уi+1). Соответствующие методы называются одношаговыми. Вычисление уi+1 в них осуществляется с использованием только одного предыдущего значения – уi, которое либо было задано по условию при i = 0, либо определилось на предыдущих итерациях. Иными словами, никаких вспомогательных процедур для реализации одношаговых методов не нужно. Поэтому их еще называют самостартующими. При k > 1 получаем многошаговые методы, более подробно поговорим о них ниже. В случае, когда функция Ф не зависит от уi+1, ее вычисление не вызывает затруднений и осуществляется по явной формуле уi+1 = уi + Ф(хi , h, уi+1–k ,… уi). Соответствующие методы называются явными. В противоположность им, методы, в которых функция Ф зависит от уi+1, называются неявными. При их реализации при каждом i возникает необходимость решения относительно уi+1 некоторого нелинейного уравнения. Пример 5.3. Требуется решить задачу Коши . Здесь f(x, y) = y2 + 1. Получим формулы некоторых простейших численных методов, используя простые (не составные) формулы численного интегрирования для аппроксимации интеграла (5.3). 1) Воспользуемся формулой левых прямоугольников, известной из курса математического анализа: = Ф(хi , h, уi). Функция Ф не зависит от уi+1 и в ней k = 1. Следовательно, получаем явный одношаговый метод уi+1 = уi + h f(хi , уi), (5.5) называемый явным методом Эйлера. Для примера 5.3 соответствующие расчетные формулы имеют вид уi+1 = уi + h [(уi)2 + 1]. 2) Воспользуемся формулой трапеций (4.5), согласно которой интеграл (5.3) имеет следующий вид: . При такой аппроксимации интеграла правая часть равенства (5.4) зависит от уi+1 и поэтому имеем неявный одношаговый метод вида уi+1 – уi = (h / 2) [ f (xi , уi ) + f (xi+1 , уi+1 ) ], (5.6) называемый неявным методом Эйлера. Согласно этому методу, в примере 5.3 для искомого значения уi+1 будем иметь квадратное уравнение, решив которое получим расчетную формулу: . В общем случае уравнение (5.6) может не иметь аналитического решения. Тогда для организации расчетов по методу (5.6) можно воспользоваться каким-нибудь численным алгоритмом решения уравнений для определения значения уi+1 при каждом i. 3) Другой способ реализации метода (5.6) без непосредственного решения уравнения относительно уi+1 состоит в соединении явного и неявного методов. Вместо точного значения уi+1 в правую часть равенства (5.6) подставляется его приближенное значение, найденное по формуле явного метода (5.5), а затем применяется формула (5.6). В результате получаем следующие расчетные формулы: = уi + h f(хi , уi); уi+1 = уi + (h / 2)[ f(xi , уi ) + f(xi+1 , ) ]. (5.7) Такая идея соединения явного и неявного методов носит название схемы прогноза и коррекции. На первом шаге по явной формуле получается величина , называемая прогнозом, которая затем корректируется с использованием формулы неявного метода. Полученное значение уi+1 считается окончательным искомым значением функции у(х) в i + 1-й точке. Метод (5.7) является одним из простейших вариантов такой схемы и называется модифицированным методом Эйлера или методом Эйлера-Коши. Подставив левую формулу (5.7) в правую, получим единую расчетную формулу, которая для примера 5.3 будет иметь вид . Далее приведены результаты расчетов перечисленных методов для примера 5.3 на отрезке [0, 1.2] c шагом h = 0.2. Для сравнения указаны точки точного решения у(х) = tg(x). Результаты расчетов примера 5.3 tg(x) – аналитич. Явный Неявный Модифицир. решение метод метод метод Графики полученных решений приведены на рис. 5.1. Видно, что наибольшая точность достигнута посредством модифицированного метода Эйлера. 2. Вычислительные формулы численных методов 4-го порядка точности Метод Рунге-Кутты. Это явный одношаговый метод. При нахождении значения уi+1 используются 4 дополнительные вычисления правой части f(x, y). , где ; ; ; . Метод Адамса-Башфорта. Это явный 4-х шаговый метод. При нахождении уi+1 используются 4 ранее найденные значения правой части fj = f(xj, yj), j = i – 3,…, i. . Метод Адамса-Моултона. Это 4-х шаговый метод типа “прогноз-коррекция”. При вычислении используется формула метода Адамса-Башфорта. Коррекция производится по формуле Моултона , , Далее приводятся результаты расчетов с помощью перечисленных методов для примера 5.3 на отрезке [0, 1.5] c шагом h = 0.15. Для сравнения указаны точки точного решения у(х) = tg(x). Лучшее совпадение с точным решением получено методом Рунге-Кутты. На рис. 5.2 приводятся графики решений. Для лучшего разрешения указаны кривые только на отрезке [0.6; 1.5] Методы 4-го порядка точности. Результаты расчетов примера 5.3 tg(x) – аналитич. Метод Метод Метод решение Рунге-Кутты Адамса-Башфорта Адамса-Моултона Правило Рунге-Ромберга. С помощью правила Рунге-Ромберга (правила двойного пересчета, аналогичного правилу Рунге для интегралов) можно подобрать значение h, обеспечивающего достижение заданной точности . Положив уh – значение искомой функции в текущей точке, вычисленное с шагом h, а yh/2 – значение, полученное при двукратном использовании тех же формул с шагом h/2, имеем следующее правило: при заданной точности  решение yh/2 следует признать удовлетворительным, если выполняется условие. , (5.8) где значение m соответствует порядку точности метода. Например, для явного метода Эйлера m = 1, для модифицированного и неявного методов Эйлера m = 2 и m = 4 для методов Рунге-Кутты и Адамса. Контрольные значения yh и yh/2 следует выбирать при таком х, для которого значение левой части (5.8) максимально. 3. Системы дифференциальных уравнений Как правило, возникающие в приложениях проблемы, приводят к необходимости решать задачу Коши не для оного дифференциального уравнения, а для систем дифференциальных уравнений вида (5.9) Здесь – искомые функции, значения которых подлежат определению при . В момент времени х = х0 задаются начальные условия , (5.10) определяющие начальное состояние физической системы, развитие которой описывается уравнениями (5.9). Введем следующие две вектор-функции , и вектор . Тогда задачу Коши (5.9), (5.10) можно записать в компактной форме , . Описанные выше численные методы решения задачи Коши для одного уравнения можно использовать и для систем уравнений первого порядка, причем форма их записи претерпевает минимальные изменения. В расчетных формулах следует лишь заменить числа yi на векторы и все функции на соответствующие вектор-функции. Например, расчетная формула метода Эйлера применительно к решению системы (5.9) принимает вид . Покоординатная запись этого соотношения выглядит так: Аналогично, формулы метода Рунге-Кутты 4-го порядка точности для систем дифференциальных уравнений первого порядка имеют вид: Решение уравнений высших порядков. Задача Коши для дифференциального уравнения n-го порядка состоит в нахождении функции y(x), удовлетворяющей при x  x0 дифференциальному уравнению (5.11) а при x = x0 – начальным условиям (5.12) Универсальным приемом решения этой задачи является ее сведение с помощью замены к задаче Коши для системы дифференциальных уравнений 1-го порядка (5.13) (5.14) Для решения задачи Коши (5.11), (5.12), приведенной к виду (5.13), (5.14), можно воспользоваться известными методами. 4. Программная реализация методов решения ОДУ 4.1. Пример программы решения системы ОДУ. Дана система из двух обыкновенных дифференциальных уравнений (5.15). Для описания принципа построения программы ее решения приведем текст программы ее численного решения модифицированным методом Эйлера. Все другие методы, в частности, 4-го порядка отличаются только количеством обращений к процедуре расчета правых частей уравнений и формулами получения результатов (прогнозируемого и окончательного). (5.15) Текст программы Uses Crt; const k = 2; { Число уравнений } type mX = array[1..k] of real; var x, x0, xN, h : real; y, f1, f2, yPR : mX; j, i, N : integer; function fAnal(x:real):real; { Аналитическое решение (чо не вставить, если есть)} begin fAnal := 2*x – 0.5*exp(x); end; procedure fun(x:real; y:mX; var f:mX); { Правые части системы } begin f[1] := –y[2] + x*x; f[2] := y[1] + exp(x); end; begin ClrScr; x0:= 0; xN:= 1; y[1]:= –0.5; y[2]:= –1.5; h:= 0.1; N:= round( (xN – x0) / h ); { Число точек решения} x := x0; writeln(x:10:5, y[1]:12:5, y[2]:10:5, fAnal(x):20:5); for i := 1 to N do begin fun(x, y, f1); for j:=1 to k do yPR[j]:= y[j] + h*f1[j]; { Прогноз } fun(x+h, yPR, f2); for j:=1 to k do y[j]:=y[j]+0.5*h*( f1[j]+f2[j]); { Результат } x := x0 + i*h; writeln(x:10:5,y[1]:12:5,y[2]:10:5,fAnal(x):20:5); readln; end; end. Таблица 5.2 Решение системы (5.14) x y1 y2 y_Aналит 0.00000 0.10000 0.20000 0.30000 0.40000 0.50000 0.60000 0.70000 0.80000 0.90000 1.00000 -0.50000 -0.35200 -0.20954 -0.07321 0.05634 0.17839 0.29215 0.39676 0.49124 0.57453 0.64545 -1.50000 -1.43724 -1.34888 -1.23432 -1.09294 -0.92401 -0.72676 -0.50032 -0.24373 0.04407 0.36425 -0.50000 -0.35259 -0.21070 -0.07493 0.05409 0.17564 0.28894 0.39312 0.48723 0.57020 0.64086 При переходе на другой метод в строке {Прогноз} записываются формулы вычисления всех прогнозируемых точек (по Адамсу-Башфорту или Рунге-Кутте), а в строках {Результат} – итоговая квадратурная формула для вычисления уi+1 с учетом полученного прогноза. Результаты расчетов приведены в таблице 5.2. 4.2. Тестовые примеры для программ решения дифференциальных уравнений 1. Решим уравнение 2-го порядка , x[ – 1, 2 ]; h = 0.375. (5.16) Данное уравнение имеет аналитическое решение у(х) = х(х2 – 1), которое является многочленом 3-й степени, а значит, методы 4-го порядка точности должны давать истинное решение. Следовательно, независимо от применяемого метода, численное решение должно в точности совпасть с аналитическим. Преобразуем уравнение (5.16) в систему уравнений 1-го порядка с помощью замены переменных: . (5.17) Начальная точка задана, следующие три точки найдем методом Рунге-Кутты 4-го порядка. Последующие точки – методами Адамса-Башфорта и Адамса-Моултона 4-го порядка (для краткости назовем их методами Адамса, т.к. для данной задачи их результаты совпадают). Функции правых частей системы (5.17) будут иметь вид: f1 (x, z1, z2) = z2; f2 (x, z1, z2) = 6 x. Таблица 5.3 Вычисление начальных точек методом Рунге-Кутты x y f1 y’ f2 k11 k12 k21 k22 k31 k32 k41 k42 –1 2 2 –6 0.75000 –2.25000 0.32813 –1.82813 0.40723 –1.82813 0.06445 –1.40625 –0.625 0.38086 0.17188 0.17188 –3.75000 0.06445 –1.40625 –0.19922 –0.98438 –0.12012 –0.98438 –0.30469 –0.56205 –0.250 0.23438 –0.81250 –0.81250 –1.5000 –0.30469 –0.56250 –0.41016 –0.14063 –0.33105 –0.14063 –0.35742 0.28125 0.125 –0.12305 –0.95313 –0.95313 0.75000 Таблица 5.4 Методы Адамса x y y’ 0.500 – 0.37500 – 0.25000 0.875 – 0.20508 1.29688 1.250 0.70313 3.68750 1.625 2.66602 6.92188 2.000 6.00000 11.00000 Результаты расчетов приведены в таблицах 5.3, 5.4. 2. Решим уравнение 2-го порядка , . (5.18) Данная задача имеет аналитическое решение у = – ln(cos(x)), которое можно использовать при тестировании программ. Преобразуем уравнение (5.18) в систему уравнений 1-го порядка: . (5.19) Таблица 5.5 Вычисление начальных точек методом Рунге-Кутты x y f1 y’ f2 k11 k12 k21 k22 k31 k32 k41 k42 0.78540 0.34657 1.00000 1.00000 2.00000 0.10472 0.20944 0.11569 0.23389 0.11697 0.23389 0.12921 0.26441 0.89012 0.46311 1.23490 1.23490 2 2.52497 0.12932 0.26441 0.14316 0.30310 0.14519 0.30310 0.16106 0.35303 0.99484 0.60763 1.53988 1.53988 3.37118 1.09956 0.78967 1.96265 1.96265 4.85184 Таблица 5.6 Метод Адамса-Башфорта x y у’ х y Y’ 1.20428 1.02434 2.58824 1.41372 1.80705 5.88182 1.30900 1.34033 3.65292 Таблица 5.7 Метод Адамса- Моултона x yпрогн у’прогн укорр у’корр уаналит 1.20428 1.02434 2.58824 1.02581 2.60813 1.02620 1.30900 1.34657 3.67281 1.35016 3.74874 1.35163 1.41372 1.83475 5.97764 1.84735 6.44233 1.85512 В таблицах 5.5 – 5.7 приведены результаты применения численных методов Рунге-Кутты, Адамса-Башфорта и Адамса-Моултона при фиксированном шаге h = /30  0.1047. С помощью правила Рунге-Ромберга подберем значение h, обеспечивающего достижение заданной точности . Для этого пересчитаем решение системы (5.18) методом Адамса-Башфорта с половинным шагом h = /60. Максимальное расхождение между yh и yh/2 получим в последней точке интервала. Имеем yh =1.807045, yh/2 =1.8450669. Отсюда  =0.002535. Полагая  = 0.0001, приходим к выводу, что точность не достигнута. Повторим дробление шага. Результаты приведены в таблице 5.8. Таблица 5.8 Применение правила Рунге- Ромберга h yh yh/2  /30 1.8070453 1.8450669 0.002535 /60 1.8450669 1.8538140 0.000583 /120 1.8538140 1.8549938 0.000079 Таким образом, требуемая точность достигается при h = /120  0.02618. Окончательное решение задачи (5.18) методом Адамса-Башфорта в сочетании с методом Рунге-Кутты приведено в следующей таблице: Таблица 5.9 Решение задачи (5.18) x у у’ Аналитическое решение – у 0.89012 0.4631143 1.23489 0.4631178 0.99484 0.6076269 1.53984 0.6079320 1.09956 0.7896605 1.96251 0.7898790 1.20428 1.0261351 2.60475 1.0261950 1.30900 1.3514021 3.73048 1.3516261 1.41372 1.8538140 6.30000 1.8551181