Анализ нелинейных резистивных схем на постоянном токе

Метод подсхем.

Большая схема разбивается на несколько многополюсных подсхем 1-го уровня, те в свою очередь на подсхемы 2-го уровня и т.д. (по мере необходимости). В результате оказывается, что вся схема и каждая подсхема состоят из небольшого числа элементов в виде многополюсников. Составляются уравнения связывающие граничные и внутренние переменные для каждого многополюсника. Анализ такой составной схемы выполняется сверху вниз. Сначала решается уравнение, составленное относительно граничных переменных схемы последней ступени, затем находятся граничные переменные подсхем предыдущей ступени и т.д. В конце концов, определяются граничные переменные подсхем первой ступени, а по уравнениям связи – все внутренние переменные. После этого расчет линейной схемы заканчивается.

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

До сих пор мы с вами рассматривали методы анализа линейных схем. Однако большинство радиотехнических устройств имеют в своем составе известные вам компоненты с нелинейной зависимостью тока от напряжения. Поэтому методы анализа линейных схем к ним в общем случае непригодны. Правда, в некоторых частных случаях, например, при работе устройства в режиме малых сигналов, нелинейные компоненты можно заменить их линейными моделями и тогда уже применять для анализа методы линейных схем. Однако и в этом случае параметры линеаризованной модели компонентов будут зависеть от режима их работы по постоянному току, т.е. от положения «рабочей точки» на их нелинейной характеристике. На постоянном токе схема становится чисто резистивной, поскольку все индуктивности становятся резисторами с нулевым сопротивлением, а емкости – резисторами с нулевой проводимостью.

4.1 Формирование узловых уравнений нелинейных схем

Рассмотрим достаточно типичный случай резистивной схемы, состоящей из линейных сопротивлений, нелинейных сопротивлений, управляемых напряжением, нелинейных источников тока, управляемых напряжением (ИТУН) и независимых источников. Это позволит нам обобщить ранее рассмотренный (см. раздел 2.1) алгоритм формирования узлового уравнения схемы на случай нелинейных цепей.

Рисунок 4.1

На рис.4.1 показана обобщенная ветвь и ее граф. Будем считать, что элемент есть либо управляемая напряжением нелинейная проводимость

, (4.1)

либо ИТУН с характеристикой

, (4.2)

где – напряжение на управляющем элементе .

Заметим, что элемент не может быть короткозамкнутой цепью, так как по условию он управляется напряжением. Уравнения (4.1) и (4.2) могут быть записаны в следующей компактной форме:

, (4.3)

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

Напомним, что закон токов Кирхгофа для узлов связной схемы может быть записан в матричном виде (3.2):

,

где A – редуцированная матрица соединений (инциденций), – вектор напряжений ветвей.

В свою очередь токи и напряжения ветвей связаны с токами i и напряжениями u на элементах следующими соотношениями (3.15-3.16):

,

,

где E и I – векторы напряжений и токов независимых источников.

Подставим выражение (4.3) в (3.16) и затем результат подставим в (3.2), откуда получим

. (4.4)

Теперь преобразование узлов (3.13):

,

подставим в (3.15) и преобразуем его к виду:

. (4.5)

Наконец, подставив напряжения на элементах (4.5) в (3.2), окончательно получим:

. (4.6)

Для схемы с N +1 узлами уравнение (4.6) является системой из N нелинейных уравнений относительно N узловых напряжений. Если обозначить вектор через x, то уравнение (4.6) может быть записано в следующей форме:

, (4.7)

где

. (4.8)

.

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

Для решения системы (4.7) обычно пользуются итерационными методами, которые состоят в том, что на каждой итерации по известному n -му приближению корня вычисляется (n+1)-е его приближение, которое лежит ближе к искомому значению точного корня . Итерационный процесс считается законченным, когда норма вектора , где – заданная погрешность вычисления. В вычислительной математике разработано большое число итерационных методов численного решения систем нелинейных алгебраических уравнений. Наиболее распространенными в схемотехническом проектировании являются метод простой итерации и модификации метода Ньютона.

4.2 Метод простой итерации для решения систем нелинейных алгебраических уравнений

В методе простой итерации исходная система уравнений (4.8) преобразуется к приведенной форме:

, (4.19)

где – итерирующая вектор-функция

, (4.20)

которую ищут в виде

. (4.21)

Подставляя (4.21) в (4.19) получаем итерационную формулу:

. (4.22)

Обычно матрица K выбирается так

. (4.23)

где – матрица Якоби системы функций () относительно переменных (x ₁, x ₂,..., x_N):

. (4.24)

При этом предполагается, что матрица неособенная.

Рисунок 4.3 – Одномерная геометрическая интерпретация неподвижной точки

Если считать некоторой точкой на оси X (как это показано на рис.4.3), то алгоритм (4.22) можно интерпретировать как средство движения точки в следующую точку . Если расстояние между последовательными парами точек уменьшается достаточно быстро, то ясно, что алгоритм будет сходится к точке , относительно которой и происходят все эти последовательные движения. С точки зрения такой интерпретации точка называется неподвижной точкой функции , а сам алгоритм иногда называют методом неподвижной точки.

При соответствующих условиях последовательность сходится к неподвижной точке:

. (4.25)

Простой критерий, который гарантирует сходимость этого алгоритма, определяется следующей теоремой.

Принцип сжатых изображений. Если есть сжатие N -мерного пространства в , т.е. имеется константа такая, что

. (4.26)

для всех , то имеет единственную неподвижную точку .

Рассмотрим подробно случай одного уравнения

. (4.27)

Его приведенная форма имеет вид

, (4.28)

а итерационный процесс строится по схеме:

, (4.29)

до выполнения условия . При этом можно рекомендовать следующую формулу для верхней оценки погрешности n -го приближения:

, (4.30)

где для , – отрезок, в котором находится корень уравнения.

В качестве одного из способов приведения уравнения общего вида (4.27) к форме (4.28) можно рекомендовать такой, который основан на следующем выборе функции :

, (4.31)

где k выбирается из условия , а знак k совпадал бы со знаком на отрезке .

Процесс решения задачи легко поддается геометрической интерпретации. Введем в плоскости декартову систему координат и построим в нем графики обеих частей уравнения – левой и правой . Абсцисса точки пересечения этих двух линий и есть решение уравнения. Геометрически процесс приближений показан на рис.4.4.

Рисунок 4.4

На (n +1)-м шаге, когда значение уже известно, вычисляется , т.е. восстанавливаем из перпендикуляр до пересечения с линией , затем из точки проводим горизонтальную прямую до пересечения с прямой , при этом абсцисса новой точки пересечения есть (n +1)-е приближение корня, т.е. .

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

и удовлетворяет условию Липшица:

, , (4.32)

где x и – любые две точки отрезка , – расстояние между ними.

Наглядный смысл соотношения (4.32) весьма простой. Функция преобразует отрезок числовой оси в отрезок , который затем без преобразования превратится в новый отрезок после операции присвоения . Отношение есть коэффициент изменения расстояния при преобразовании оператором j. Величина p ограничена числом L. Условие говорит о том, что при преобразовании происходит уменьшение отрезка , следовательно операция является сжатием. Если обозначить погрешность вычисления корня на n -м шаге итерации символом , то скорость сходимости можно определить приближенным соотношением:

. (4.33)

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

4.3 Метод Ньютона

Согласно методу Ньютона последовательные приближения для корней системы нелинейных уравнений находятся по формуле:

, (4.34)

где – поправки вектора, – матрица Якоби системы функций () относительно переменных (x ₁, x ₂,..., x _N).

За нулевое приближение можно найти грубое приближение искомого корня.

В случае одного уравнения с одним неизвестным матрица Якоби вырождается в обычную производную и выражение (4.34) преобразуется к виду:

, (). (4.35)

Правило (4.35) имеет простой геометрический смысл. В плоскости с системой координат xOy построим график функции (рис.4.5).

Рисунок 4.5

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

Можно показать, что метод Ньютона имеет квадратичную сходимость:

, где . (4.36)

Из равенства (4.36) следует, что при будет иметь место неравенство , т.е. последующее приближение будет ближе к искомому корню, чем предыдущее, а следовательно итерационный процесс сходится. Таким образом, при правильном выборе начального приближения метод Ньютона всегда сходится.

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

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

4.4. Решение узловых уравнений методом Ньютона

Применим n -мерный метод Ньютона для решения нелинейного узлового уравнения:

. (4.37)

Матрица Якоби от имеет вид:

. (4.38)

Следовательно, итерационный процесс Ньютона строится по такой формуле:

. (4.39)

Установим физический смысл полученного уравнения, для чего перепишем его в следующем виде:

. (4.40)

и сравним с узловым уравнением для линейной схемы

, (3.18)

в котором

, .

Уравнения (4.40) и (3.18) становятся очень похожими, если ввести следующие обозначения

; (4.41)

; (4.42)

; (4.43)

. (4.44)

В таком случае (4.40) принимает вид:

,

где и – матрицы дифференциальных проводимостей ветвей и узлов, а и – приращение узловых источников токов и узловых напряжений на (n+1)-й итерации, которые показаны графически на рис.4.6 для k -й ветви.

Рисунок 4.6

Отметим также, что в частном случае, когда все сопротивления линейны, сводится к обычной матрице проводимостей ветвей , и поэтому не зависит от номера итерации i.

5 МОДЕЛИРОВАНИЕ электрических схем
ВО ВРЕМЕННОЙ ОБЛАСТИ

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

5.1. Математическая модель схемы для временной области.

Метод переменных состояния

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

Рисунок 5.1

Записывая ЗТК для 1-го и 2-го узла, получаем систему интегро-дифференциальных уравнений следующего вида:

;

или в преобразованном виде:

;

.

Решение интегро-дифференциальных уравнений есть достаточно сложная математическая задача, решение которой наталкивается на ряд вычислительных трудностей.

Значительно проще решать системы уравнений, которые не содержат интегральных слагаемых, а имеющие в своем составе только производные неизвестных функций, т.е. дифференциальные уравнения. Методы решения дифференциальных уравнений достаточно хорошо разработаны, в том числе и численные методы, ориентированные на применение ЭВМ. Поэтому в последние годы получил широкое распространение метод составления уравнений цепи, который приводит к модели схемы в виде системы обыкновенных дифференциальных уравнений. В этом методе, который получил название метода переменных состояния, неизвестными в системе уравнений являются так называемые переменные состояния, а именно напряжения на емкостях и токи в индуктивностях анализируемой электрической цепи. Эти переменные определяют как общий запас энергии в цепи, так и энергию, запасенную в каждом реактивном элементе цепи, т.е. определяют электромагнитное состояние цепи. Если значения этих неизвестных, т.е. переменных состояния, будут найдены, то после этого легко находятся токи и напряжения на зажимах любых других элементов цепи.

Действительно, если в некоторый момент времени найденное напряжение на p- й емкостной цепи равно , то по теореме замещения этот элемент можно заместить источником напряжения с тем же задающим напряжением. Аналогично для того же момента времени q -ю индуктивность можно заменить источником тока с задающим током , значение которого известно. Если подобную замену выполнить для всех реактивных элементов цепи, то в результате образуется резистивная электрическая цепь, токи и напряжения в которой могут быть найдены методами анализа резистивных электрических цепей.

Для примера рассмотрим ту же цепь (рис.5.1), которую для удобства приведем к виду, показанному на рис.5.2, а. Если, например, найдены переменные состояния и в этой цепи, то после указанной замены могут быть легко рассчитаны токи и напряжения в резистивной цепи, схема которой изображена на рис.5.2, б. Что касается определения самих переменных состояния, то они могут быть определены из уравнений, которые составляются на основе применения законов Кирхгофа.

а) б)

Рисунок 5.2

Действительно, в соответствии с 1-м законом Кирхгофа

,

поэтому

.

В соответствие же со 2-м законом Кирхгофа

.

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

. (5.1)

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

;

; (5.2)

................................

.

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

, .

Матричная форма записи уравнений (5.2) имеет вид:

, или , (5.3)

где – матрица параметров цепи; – матрица-столбец переменных состояния; ; B – квадратная матрица; – матрица-столбец воздействий в виде токов и напряжений независимых источников.

Для рассматриваемого примера:

; ; ; . (5.4)

Таким образом, выбор переменных состояниям в виде токов в индуктивностях и напряжений на емкостях цепи в качестве неизвестных, подлежащих определению, приводит к тому, что:

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

б) сами переменные состояния характеризуют энергию, запасенную в реактивных элементах анализируемой цепи.

5.2. Решение уравнений переменных состояния линейных цепей

во временной области

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

с начальным условием ,

или в более общем виде:

, (5.16)

где , – матрица-столбец неизвестных (переменных состояния), – столбец воздействий, A и B – квадратные матрицы.

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

. (5.17)

Если , то говорят, что линейная система является правильной, т.е. для нее

.

Решение уравнения (5.3) будем искать в виде:

. (5.18)

Неизвестной здесь является функция . Найдем ее, для чего подставим (5.18) в уравнение (5.3):

,

откуда следует:

. (5.19)

Интегрируем последнее равенство от до t:

(5.20)

получаем:

. (5.21)

Из уравнений (5.18) и (5.21) имеем

, (5.22)

Найдем . Устремляя в (5.22), имеем:

.

откуда

. (5.23)

Окончательно:

. (5.24)

. (5.25)

Равенство (5.24) является точным решением дифференциального уравнения (5.16), однако его не очень удобно использовать в прямую для вычисления на ЭВМ. Обычно вычисляют для фиксированных значений , где – целые числа,– выбранный временной интервал. Для этого равенство (9) сводят к так называемому разностному уравнению. Рассмотрим один из способов такого сведения.

Пусть функция возбуждения задана в виде, изображенном на рис.5.5.

Рисунок 5.5.

Рисунок 5.5

Разбиваем весь интервал интегрирования на равные интервалы D, тогда точка начала каждого интервала определится, как . Значение искомой функции в конце интервала можно определить из (5.24):

. (5.26)

Для облегчения операции интегрирования в (5.26) аппроксимируем функцию кусочно-постоянной функцией , которая в пределах каждого интервала будет постоянной и равной значению в начале этого интервала (см. рис.5.5), т.е.:

.

Подставляя в (5.26) и используя тождество

,

имеем

. (5.27)

Подставляя (5.27) в (5.26), находим

, (5.28)

которое и является искомым разностным уравнением.

Полученное равенство по существу является формулой, легко поддающейся вычислению на ЭВМ. Входящие в него матричные выражения легко вычисляются с помощью разложения в ряд Тейлора:

5.3. Модель электрической схемы с нелинейными элементами

Рассмотрим электрическую цепь, состоящую из линейных и нелинейных резисторов, емкостей и индуктивностей, а также независимых источников тока и напряжения. Будем полагать, что нелинейные емкости управляются напряжением, а нелинейные индуктивности – протекающими через них токами. Кроме того, цепь имеет такую же топологию, которое допускает построение правильного дерева, где все емкости и управляемые напряжением резисторы относятся к ветвям дерева, а все индуктивности и управляемые током проводимости – к дополнению дерева.

При таких условиях модель схемы будет представлять собой систему нелинейных ОДУ, которую можно записать в матричном виде:

, (5.37)

где

. (5.38)

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

, (5.39)

для решения которого надо задать начальные условия , и найти неизвестные матрицы из уравнения (5.37).

Для решения (5.37) воспользуемся методом Ньютона. В (5.37) обозначим:

; ,

тогда уравнение (5.37) приводится к виду:

,

для решения которого можно организовать процесс Ньютона:

который для неизвестных , записывается в следующем виде:

. (5.40)

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

Шаг 1. .

Начальные условия: , .

1. Находим из (5.38).

2. Находим значения , из решения нелинейного алгебраического уравнения (5.37) по формуле (5.40)

3. Вычислим V из (5.35)

4. Решаем дифференциальное уравнение и находим значения и в точке

Шаг 2. .

, .

Все действия в шагах (1-4) повторяются и находятся и для момента .

Затем следуют шаги 3 … N, пока не будут определены все значения функции на интервале Т.

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

Уравнение состояния в компактной нормальной форме имеет вид:

. (5.41)

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

1) Разложение искомой функции в ряд Тейлора;

2) Использование полиномиальной аппроксимации.

Во всех методах значения функции ищутся в дискретных точках . Для упрощения положим, что все эти точки расположены на всем интервале равномерно с шагом h. При таком допущении интервал времени определяется соотношением:

, где < T.

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

5.6. Методы Эйлера для численного решения обыкновенных дифференциальных уравнений

Методы Эйлера относится к одношаговым методам, использующим разложение искомой функции в ряд Тейлора. Рассмотрим их на примере решения одного обыкновенного дифференциального уравнения вида:

.

Требуется найти функцию на интервале [] с начальным условием

, .

Выберем достаточно малый шаг h по аргументу t и построим на [] систему равноотстоящих точек , ().

Найдем значение функции в точке , для чего представим ее в виде разложения в ряд Тейлора в окрестности точки

.

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

.

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

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

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

носит название формулы Эйлера.

Геометрическая интерпретация итерационного процесса вычисления с помощью формулы Эйлера может быть выяснена с помощью рис.5.7.

Рисунок 5.7.

В начальной точке нам известно точное значение функции . На первом шаге мы заменяем искомую функцию отрезком прямой, касательной к в начальной точке отрезка []. Найденная величина отличается от точного значения функции на некоторую величину .

Следующий шаг выполняется на отрезке аргумента [] из точки () в точку () также по прямой линии, параллельной касательной к в точке , и т.д. Таким образом, искомая функция заменяется кусочно-линейной ломаной кривой с изломами в точках (), (),... Нетрудно увидеть, что на каждом шаге в качестве начальной точки используется приближенное значение функции в предыдущей точке, что приводит к накоплению ошибки по мере движения вдоль отрезка [].

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

, (),

, (),

то имеет место следующая оценка погрешности

.

Последняя формула говорит о том, что погрешность вычислений пропорциональна величине шага h, поэтому действенным способом повышения точности является уменьшение шага.

Несомненным достоинством метода Эйлера является его простота, а недостатками – малая точность и накопление ошибок. Существуют две модификации метода Эйлера, которые в некоторой степени уменьшают эти недостатки, повышая точность вычисления.

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

,

где

,

.

Объединяя последние равенства, окончательно имеем

.

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

,

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

и подставить ее в правую часть

,

Объединяя равенства последние равенства, окончательно находим формулу, представляющую второй улучшенный метод Эйлера:

.

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

5.6.2. Методы Тейлора (общий случай)

В методах Эйлера для построения вычислительной схемы использовалась только первая производная искомой функции. Рассмотрим более общий случай, когда при разложении искомой функции используется произвольное число членов ряда Тейлора. Численное решение ОДУ методом рядов Тейлора применительно к одному дифференциальному уравнению, приведенного к каноническому виду:

при начальном условии .

Разложим неизвестную функцию в ряд Тейлора в окрестности точки , полагая величину шага h малой:

. (5.42)

Сделаем в уравнении (5.42) следующие замены:

и , ,…., , …

и усечем ряд:

Введем обозначение:

(5.43)

тогда

. (5.44)

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

Обозначим погрешность метода величиной e:

В математике доказано, что эта погрешность зависит от величины шага и порядка метода :

.

Метод Тейлора первого порядка, по сути, является явным методом Эйлера (или просто методом Эйлера), в чем нетрудно убедиться, в (5.43) и (5.44) :

.

В методе Тейлора второго порядкав (5.44) следует положить , тогда:

,

, (5.45)

где ; .

Аналогично можно получить формулы для метода Тейлора третьего, четвертого и т.д. порядков. Очевидный недостаток всех этих методов при >1 – необходимость вычисления частных производных, что является трудной задачей и является источником ошибок, особенно когда не выражается в аналитической форме.

5.6.3. Методы Рунге – Кутта

Математики Рунге и Кутт в начале ХХ века нашли метод, позволяющий исключить необходимость вычисления частных производных и имеющих тот же порядок точности, что и алгоритм Тейлора.

Основная идея их метода состоит в том, что функция заменяется другой функцией не содержащей частных производных от , так что:

, (5.46)

где некая константа, не зависящая от .

С учетом (5.46) модифицированный метод

(5.47)

имеет тот же порядок методической ошибки, что и метод Тейлора. Этот метод получил название метода Рунге – Кутта порядка .

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