Лабораторная работа 6 Численное дифференцирование

Работа добавлена на сайт samzan.net:

PAGE  15

НГТУ, факультет РЭФ, кафедра ППиМЭ

Лабораторный практикум «Информатика-3» Алгоритмы и программы вычислительных задач микро- и наноэлектроники

Лабораторная работа № 6

Численное дифференцирование. Численное интегрирование обыкновенных дифференциальных уравнений

© Н.В.Усольцев       Редакция 2007 г.       Вариант № 1 от 16.11.08

1. Введение

При решении различных вычислительных задач, часто возникает необходимость расчета не только интегралов, но и значений производных (первой, второй и более старших порядков), различных функций, том числе, и заданных алгоритмически или таблично. Вычисление производных основывается на их представлении в виде конечных разностей.

На этом же основывается и численное интегрирование дифференциальных уравнений. Обыкновенные дифференциальные уравнения (ОДУ) содержат под знаком дифференцирования неизвестную функцию одной переменной , которая и является искомым решением. В зависимости от порядка старшей производной ОДУ могут быть первого, второго, третьего и более высоких порядков. Несколько ОДУ могут образовывать систему одновременно интегрируемых уравнений. В этом случае результатом решения системы является совокупность функций .

Численное интегрирование ДУ отличается от аналитического, прежде всего тем, что результат получается в табличном численном виде. Кроме того, они обладают большей универсальностью и относительной простотой.

В данной работе изучаются методы численного дифференцирования, а также численные методы интегрирования ОДУ первого и второго порядков и систем ОДУ 1-го порядка.

2. Методы численного дифференцирования

Методы численного дифференцирования функции одной переменной  основаны на формальном определении производной:

Самыми простыми формулами для вычисления первой производной  являются двухточечные формулы дифференцирования вперед и назад:

      (1)

      (2)

Они дают погрешность примерно одинаковой величины, но разного знака и поэтому обычно используют среднее значение производных, получаемых по этим формулам:

  (3)

Эта формула называется трехточечной, т.к. в ней используется значение функции в трех точках: ,  и , причем последнее в ней присутствует неявно.

Определение второй производной как производной от первой

приводит к формуле

     (4)

При использовании формул численного дифференцирования результат зависит от правильного выбора дифференцирующего приращения аргумента .  Оно должно быть достаточно малым.

Приведенные выше формулы предполагают возможность расчета производной при возможности вычисления функции в любой точке x. Если же функция задана таблично и шаг аргумента между точками достаточно мал, для вычисления производных в -той точке можно воспользоваться аналогичными (1 - 4) формулами

- формула дифференцирования вперед;   (5)

- формула  дифференцирования назад;   (6)

- трехточечная формула.     (7)

3. Методы численного интегрирования ОДУ 1-го порядка

ОДУ 1-го порядка в большинстве случаев оно допускает явное выражение производной (форма Коши)

       (8а)

Для нахождения его однозначного решения необходимо начальное условие

          (8б)

Дифференциальное уравнение (8а) и условие (8б) составляют задачу Коши, которая имеет однозначное решение.

В дальнейшем полагается, что решение необходимо найти на интервале , хотя можно находить его и в сторону уменьшения .

Численные методы интегрирования основаны на том, что интервал поиска решения, начиная от точки , разбивается на элементарные участки длиной , которая называется длиной шага или просто шагом. Будем использовать ее величину постоянной. Длина шага должна быть достаточно мала; чем она меньше, тем точнее будет результат. Значения аргумента в точках разбиения образуют последовательность:

.

В каждой из этих точек решение задачи Коши  имеет значения , которые пока неизвестны. Их и нужно найти. Как это делается, будет показано ниже, а пока следует отметить, что результат численного интегрирования будет представлять собой функцию  в табличной форме. Далее она может быть приведена к аналитической форме методами аппроксимации, но в большинстве случаев можно ограничиться и табличным решением.

После разбиения интервала интегрирования, можно записать ОДУ в k-той точке в виде

      (9)

а затем представить в конечно-разностной форме.

Явный метод Эйлера

Если подставить формулу дифференцирования вперед (5)  в (9), а затем выразить , то получим формулу численного интегрирования явного метода Эйлера:

        (10а)

где  k = 0, 1, 2, 3, …- индекс шага.

Она используется вместе с очевидной формулой

         (10б)

Процесс численного интегрирования явным методом Эйлера представляет собой расчет таблицы значений  и  по формулам (10), начиная с  и  и с использованием  и  в качестве стартовых значений.

Численное интегрирование обычно заканчивается, когда пройден заранее заданный диапазон  или значения функции перестали изменяться.

Неявный метод Эйлера

Если в представлении производной  использовать формулу дифференцирования назад (6), то получается формула численного интегрирования неявного метода Эйлера:

;      (11a)

       (11b)

Она представляет собой алгебраическое уравнение относительно искомой величины . При общем нелинейном виде правой части ОДУ , оно не допускает аналитического решения и может решаться численными методами (например, простых итераций или Ньютона). Иными словами, неявный метод Эйлера предполагает на каждом шаге интегрирования решение нелинейного уравнения и, поэтому, он более трудоемкий, чем явный. Однако существуют причины его применения, которые будут разъяснены далее.

Результат интегрирования также представляет собой функцию в табличном виде.

На точность представления производной в конечно-разностном виде (5) или (6) влияет величина шага . С ее уменьшением погрешность интегрирования также уменьшается, но до определенных пределов. При очень малых шагах начинает играть роль погрешность вычислений и накопление суммарной погрешности, т.к. совершается большое число шагов.

При выборе величины шага  ее следует соизмерять с величиной        ,         (12)

которая называется постоянной интегрирования. Это название неточно, т.к. при нелинейной правой части ОДУ  она не постоянна и может изменяться в диапазоне интегрирования весьма существенно.

Для большого класса ОДУ величина шага

         (13)

обеспечивает погрешность не более 1…2%. Обычно такую величину шага и считают оптимальной.

Метод Рунге-Кутты

Явный и неявный методы Эйлера используют в представлении производной только линейное соотношение между приращением аргумента и приращением функции. Если же представить

,

то получаются следующие формулы численного интегрирования метода Рунге-Кутты:

 ;  ;

где   ;

 ;        (14)

 ;

 ;

 .

Метод Рунге-Кутты обеспечивает очень высокую точность: даже при длине шага  на гладких функциях погрешность не превышает сотых долей процента.

Метод относится к явным, т.к. на каждом шаге не нужно решать уравнений; необходимо выполнить последовательные вычисления по формулам (14).

Устойчивость численного интегрирования

Процесс численного интегрирования может быть неустойчивым. Числовая неустойчивость представляет собой катастрофическое накопление ошибки, при котором результат интегрирования с каждым шагом все сильнее отклоняется от истинного решения. Неустойчивость проявляют только явные методы при больших шагах интегрирования.

Для явного метода Эйлера условием устойчивости является

          

Неявный метод устойчив при любой величине шага; в это заключается его преимущество, которое особенно используется при интегрировании систем ОДУ 1-го порядка.

4. Методы численного интегрирования систем ОДУ 1-го порядка

Система ОДУ 1-го порядка содержит  уравнений и решением его является набор из  функций одной переменной:

 .

Линейная относительно производных система ОДУ 1-го порядка может быть записана в виде:

 (15)

В общем случае коэффициенты  также могут зависеть от .

Если ввести векторные обозначения  и , то система (13) может быть представлена в векторно-матричном виде

        (16)

Если определитель матрицы  не обращается в нуль, то умножая обе части (16) на матрицу , получаем явный вид системы относительно вектора производных:

          (17а)

или, в развернутом виде:

          

Для однозначного решения  система должна иметь набор начальных условий:

 

или в векторной форме:

          (17б)

Аналогия векторной записи системы ОДУ (17а) и отдельного уравнения (8а) позволяет легко записать формулы численного интегрирования

•  явного метода Эйлера:

     (18)

•  неявного метода Эйлера:

     (19)

•  метода Рунге-Кутты:

 ;  

где   ;

 ;       (20)

 ;

 ;

 .

Входящие в формулы (18) и (20) векторные операции легко выполняются. Выражение (19) представляет собой систему нелинейных алгебраических уравнений относительно вектора . Она может решаться методом простых итераций или Ньютона.

Выбор величины шага для численного интегрирования системы осуществляется на основе анализа набора постоянных интегрирования . Они находятся следующим образом. Сначала для системы рассчитывается матрица Якоби, элементы которой есть производные правых частей уравнений по неизвестным функциям:

          (21)

Затем необходимо найти собственные значения матрицы Якоби, т.е. корни характеристического уравнения

 ,        (22)

где   - символ определителя,

 - единичная матрица,

 - неизвестная, относительно которой осуществляется решение.

Левая часть характеристического уравнения (22) представляет собой полином -ой степени , поэтому имеет  корней, в общем случае, комплексных: ,      , где  - мнимая единица.

Постоянные интегрирования  определяются по собственным значениям:

        (23)

При использовании явного метода Эйлера величина шага  должна удовлетворять условию устойчивости

 ,         

а чтобы обеспечить и достаточную точность интегрирования

 ,        (24)

где  - минимальная по абсолютной величине .

Следует иметь в виду, что , в строгом смысле, не являются постоянными. Если правые части системы  нелинейные по , то матрица Якоби, ее собственные значения  и “постоянные интегрирования”  изменяются на каждом шаге, причем, иногда на порядки. Все это порождает сложности при интегрировании систем.

Системы, постоянные интегрирования которых имеют значительный разброс, называются жесткими. При интегрировании жестких систем обычно используется неявный метод Эйлера, как устойчивый при любой величине шага.

5. Методы численного интегрирования ОДУ 2-го порядка

В большинстве физических и технических задач встречаются линейные ОДУ 2-го порядка, имеющие общий вид:

     (25)

Функция , ее первая и вторая производные входят в него в линейной комбинации; коэффициенты , ,  и правая часть уравнения  могут быть любыми функциями от .

Для решения уравнения (25) необходимо два граничных условия. В зависимости от их вида различают

•  задачу Коши:

;  ;

•  краевую задачу:

;  ;

где  и  - границы интервала поиска решения.

Для численного интегрирования краевой задачи интервал поиска решения  разбивается на  элементарных интервалов с шагом :

 ;         

Получается набор  узловых точек

;  ;

Искомая функция  имеет значения (пока неизвестные во всех точках, кроме крайних)   

;

Уравнение (25) записывается в произвольной -той узловой точке

 ,     (26)

где  ; ; ; .

Подставив в (26) представление производных в виде конечных разностей

;

 

и перегруппировав члены, получаем:

 ,      (27)

где ; ; ;

При разворачивании индекса  от 1 до  выражение (27) проявляет свою математическую сущность – это система линейных алгебраических уравнений относительно :

 .

Перенося в правые части известные члены  и , видим, что матрица системы имеет трехдиагональный вид:

 

В сокращенном векторно-матричном виде

         (28)

Вектор  представляет собой вектор граничных условий. Он имеет две ненулевых компоненты: первую и последнюю.

Однократное решение системы (28) позволяет получить искомые значения функции  в узловых точках. Решение проводится методом прогонки.

Некоторые дифференциальные уравнения, используемые в работе

Для проверки методов численного интегрирования удобно использовать дифференциальные уравнения, решение которых можно получить аналитически и сравнить его с численным. В качестве такого простейшего ОДУ 1-го порядка может быть взято уравнение

;  .       (29)

Его решение имеет вид:  

        (30)

Таким уравнением описывается, например, процесс зарядки  емкости C через резистор R при ступенчатом изменении входного напряжения от 0 до E (рис.1). Уравнение, полученное на основе закона токов Кирхгофа, имеет вид:

 или

,  где .      

Его решение имеет вид

      

Для интегрирования систем ОДУ 1-го порядка в задании к работе предлагается система, описывающая переходные процессы в схеме (рис.2) с двумя внутренними узлами, что порождает два уравнения. Уравнения, составленные для нее по закону токов Кирхгофа, имеют вид:

Она легко приводится к виду (14)

;  (31)

или

 

Умножая обратную матрицу  на левую и правую части уравнения, получим:

,      (32)

где  - матрица Якоби системы, а .

При выборе шага интегрирования необходимо найти собственные значения матрицы Якоби и постоянные интегрирования (постоянные времени). Характеристическое уравнение (22) будет представлять собой квадратное уравнение:

;

или

,       (33)

где   - определитель матрицы Якоби.

Это уравнение легко решается аналитически и позволяет найти собственные значения и постоянные времени схемы для выбора шага.

8. Программное обеспечение к работе

Для вычисления первой и второй производных можно использовать функции модуля dif.cpp. Они реализуют формулы (1) – (4).

Программы для численного интегрирования ОДУ и систем ОДУ 1-го порядка объединены в модуль DifEqu1 (файл difequ1.cpp).

Для интегрирования отдельных уравнений используется функция IntegroOne(). Она реализует явный или неявный методы Эйлера, а также метод Рунге-Кутты (в зависимости от входной переменной – ключа метода).  Правая часть ОДУ должна быть оформлена в виде функции с заголовком double FuncName(double y, double x). Она обычно располагается в одном файле с управляющей программой.

Функция IntegroOne()обращается к функциям, выполняющим один шаг интегрирования различными методами:

•  EilerExpOne() – реализует явный (explicit) метод Эйлера;
•  EilerImpOne() – реализует неявный (implicit) метод Эйлера;
•  RungeKutOne() – реализует метод Рунге-Кутты.

Функция EilerImpOne() имеет дополнительный параметр – ключ key, задающий метод решения уравнения на каждом шаге интегрирования. Могут быть использованы метод простых итераций или метод Ньютона. Они “встроены” в функцию, т.е. не требуют подключения модуля Root. Метод Ньютона следует использовать при больших шагах интегрирования, когда метод простых итераций может расходиться. Признаком расходимости являются нулевые значения Yk, возвращаемые EilerImpOne().

Для интегрирования систем ОДУ 1-го порядка используется функция IntegroSys(). Она устроена аналогично IntegroOne(). Правая часть системы ОДУ должна быть оформлена в виде функции с заголовком  void FuncName(double x, Vector Y, Vector F).

Функция IntegroSys()также обращается к функции, выполняющей один шаг интегрирования явным методом Эйлера EilerExpSys() .

Программы для численного интегрирования ОДУ 2-го порядка объединены в модуль DifEqu2 (файл difequ2.cpp).

Для численного интегрирования краевой задачи с линейным ОДУ 2-го порядка вида (25) используется функция de2abLin() из модуля Difequ2. Конкретное уравнение задается функциями его коэффициентов ,  ,  и правой части . Эти функции должны иметь заголовки типа double FuncName(double x) и обычно располагаются в одном файле с управляющей программой.

В модуль DifEqu2 входит также функция LineProgon(), реализующая метод прогонки.

Списки входных и выходных параметров функций предлагается изучить по их входным комментариям самостоятельно.

9. Задание к работе

1.  Для функции  в точке x = 0 выполнить расчет относительной погрешности вычисления первой производной по формуле дифференцирования вперед (1) и трехточечной формуле (3), а также второй производной по формуле (4)  при различных значениях дифференцирующего приращения dx. Величину приращения логарифмически изменять от 10-10 до единицы. Расчеты представить на графике в логарифмических масштабах по осям. Объяснить полученные результаты.

1.  Проанализировать функции модуля difequ1.cpp и убедиться в их соответствии с описанными выше алгоритмами.

1.  Выполнить численное интегрирование  задачи Коши:  ;   на интервале (0…3) с шагом
   •  явным методом Эйлера;
   •  неявным методом Эйлера;
   •  модифицированным методом Эйлера;
   •  методом Рунге-Кутты.

Построить все решения на одном графике вместе с точным решением (файл exp3.dat). Выполнить сравнение результатов и сделать вывод о погрешности различных методов.

1.  Выполнить численное интегрирование этой же задачи явным методом Эйлера с шагом . Представить на одном графике это решение, решение, полученное в предыдущем пункте с шагом  и точное аналитическое решение. Сделать вывод о зависимости точности решения от величины шага.

1.  Выполнить численное интегрирование явным методом Эйлера этой же задачи на интервале (0…15) с шагами:

•  
•  
•  
•  ,

а также неявным методом Эйлера с шагом . Представить все решения на одном графике. Туда же поместить точное аналитическое решение из файла exp15.dat. Проанализировать результат и сделать выводы об устойчивости численного интегрирования.

1.  Выполнить численное интегрирование явным методом Эйлера системы (35), описывающей процессы в схеме на рис. 2, при заданных в таблице значениях параметров (по вариантам). Для правильного выбора шага интегрирования по формуле (23) необходимо предварительно рассчитать постоянные времени по формулам (38) и (20). Расчет можно выполнить с помощью таблицы, находящейся в файле tau.xls. При расчетах не нужно переводить сопротивления в Омы, а емкости – в Фарады, т.к. 1кОм1пФ = 1 нс. Основное изменение потенциалов при заданных параметрах будет происходить в диапазоне единиц, десятков или сотен наносекунд.

№ варианта	Е, В	R1, кОм	R2, кОм	С1, пФ	С2, пФ	С3, пФ
1	10	1	2	1	2	0.5
2	5	5	10	5	2	1
3	25	6	3	10	15	5
4	1	10	4	2	4	1
5	15	2	1	4	5	2
6	5	10	15	5	10	2
7	50	3	5	10	5	5
8	20	5	10	4	2	1

1.  Проанализировать функции модуля difequ2.cpp и убедиться в их соответствии с описанными выше алгоритмами.

1.  Выполнить решение краевой задачи с линейным ОДУ 2-го порядка. Решения для разных значений параметра ν построить на одном графике.
   1.  ;   y(-1) = -1; y(1) = 1;  

для ν = 1, 3, 5, 7;

1.  ;    y(0) = -1; y(2) = 1;

для ν = 2, 4, 6, 8;

1.  ;  y(0)=1; y(10)=-0.2459;

для ν = 0, 1, 2, 3;

1.  ;   y(0)=1; y(5)=1;

для ν = 2, 3, 4, 5;

1.  ;  y(-1) = 1; y(1) = 1;

для ν = 2, 4, 6, 8;

1.  ; y(0) = 0; y(1) = 1;

для ν = 1, 2, 3, 4;

1.  ; y(0)=0; y(10)=1;

для ν = 1, 2, 3, 4;

1.  ;  ;

для ν = 3, 4, 5, 6;

1. а третьей группы из малообеспеченных семей; студентыучастники боевых действий; студенческие семьи г
2. Гомельский государственный университет имени Франциска Скорины Исторический факультет Кафедра
3. Отражение результатов в бухгалтерском учете 6 Учет товаров в комиссионных магазинах 7 Учет реализации то
4. все ответы верны 2
5. Тема 5- Школа научного управления 1885 1920 гг
6. Термохемилюминесцентный иммуноанализ
7. кишечного тракта способно значительно облегчить протекание голодания особенно у начинающих
8. Повесть временных лет относит призвание варягов норманнского князя Рюрика 882г
9. неформальные частные предположения о природе событий предметов и ситуаций с которыми мы сталкиваемся в жи
10. Бисквитнокремовый Норма расхода кг-т
11. практическая конференция ИННОВАЦИОННОЕ РАЗВИТИЕ СОВРЕМЕННОЙ НАУКИ 31 ЯНВАРЯ 2014Г
12. Менталитет русской культуры
13. а актілер аны~тамалар ма~л~мат хаттар мен т~сінік хаттар арыздар т
14. і Довгий шлях розвитку пройшло людство від первісного вогнища в давні часи і до сучасного місця за обідні
15. Размножение водорослей
16. Психодраматический подход к проблемам человека
17. на тему- ldquo;КАТЕГОРИИ ДИАЛЕКТИКИrdquo; МОСКВА 1998 Н
18. . Микроэкономика ~ представляет собой раздел экономической теории изучающей поведение потребителей и фирмы.
19. Развитие транспорта в Москве
20. тематического аппарата квантовой механики

Материалы собраны группой SamZan и находятся в свободном доступе