Вариант 1

1.1 1.2

1.3 1.4 1.5

Вариант 2

2.1 2.2

2.3 2.4 2.5

Вариант 3

3.1 3.2

3.3 3.4 3.5

Вариант 4

4.1 4.2

4.3 4.4 4.5

Вариант 5

5.1 5.2

5.3 5.4 5. 5

Вариант 6

6.1 6.2

6.3 6.4 6.5

Вариант 7

7.1 7.2

7.3 7.4 7.5

Вариант 8

8.1 8.2

8.3 8.4 8.5

Вариант 9

9.1 9.2

9.3 9.4 9.5

Вариант 10

10.1 10.2

10.3 10.4 10.5

Вариант 11

11.1 11.2

11.3 11.4 11.5

Вариант 12

12.1 12.2

12.3 12.4 12.5

Вариант 13

13.1 13.2

13.3 13.4 13.5

Вариант 14

14.1 14.2

14.3 14.4 14.5

Вариант 15

15.1 15.2

15.3 15.4 15.5

Вариант 16

16.1 16.2

16.3 16.4 16.5

Вариант 17

17.1 17.2

17.3 17.4 17.5

Вариант 18

18.1 18.2

18.3 18.4 18.5

Вариант 19

19.1 19.2

19.3 19.4 19.5

Вариант 20

20.1 20.2

20.3 20.4 20.5

Вариант 21

21.1 21.2

21.3 21.4 21.5

Вариант 22

22.1 22.2

22.3 22.4 22.5

Вариант 23

23.1 23.2

23.3 23.4 23.5

Вариант 24

24.1 24.2

24.3 24.4 24.5

Вариант 25

25.1 25.2

25.3 25.4 25.5

Вариант 26

26.1 26.2

26.3 26.4 26.5

Вариант 27

27.1 27.2

27.3 27.4 27.5

Вариант 28

28.1 28.2

28.3 28.4 28.5

Вариант 29

29.1 29. 2

29.3 29.4 29.5

Вариант 30

30.1 30.2

30.3 30.4 30.5 .

ТЕМА. Системы нелинейных уравнений

При решении задач моделирования поведения химических систем достаточно часто приходится решать системы уравнений, нелинейных по отношению к переменным. Системыn линейных уравнений с n неизвестными x1, x2, ..., xn  в общем случае принято записывать следующим образом:

где F1, F2,…, Fn – любые функции независимых переменных, в том числе и нелинейные относительно неизвестных.

Как и в случае систем линейных уравнений, решением системы является такой вектор (или векторы) (X*), который  при подстановке обращает одновременно все уравнения системы в тождества.

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

Рассмотрим несколько простейших итерационных методов решения систем нелинейных уравнений, а именно, метод простой итерации, метод Зейделя и метод Ньютона.

Метод простой итерации

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

Выбирая затем вектор начального приближения

,  

подставляют его в преобразованную систему уравнений. Из первого уравнения получают новое приближение к первой переменной, из второго – второй и т. д. Полученное уточненное значение переменных снова подставляют в эти уравнения и т.д.Таким образом, на (i+1)-м шаге итерационной процедуры имеем

Метод Зейделя

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

Метод Ньютона-Рафсона

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

Рассмотрим метод на примере системы двух уравнений с двумя неизвестными:

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

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

после пренебрежения всеми членами, кроме линейного:

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

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

Запишем для функции F1 2-х переменных линейную часть разложения в ряд Тейлора в окрестности выбранной точки

для второго уравнения, аналогично

Если значения переменных x1 и x2 являются решением, то оба уравнения системы должны обратиться в ноль, поэтому полученные разложения приравниваем нулю.

Для краткости записи введем следующие обозначения:

 - приращение i-ой переменной  

- значение первой частной производной функции Fj по переменной xi при значении переменных 

 – значение  j-ой функции при соответствующих значениях переменных, то есть невязка j‑го уравнения. 

Получим систему линейных уравнений 2 x 2 относительно приращения переменных

Или, в матричной форме,

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

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

Итерационная процедура далее продолжается аналогично.

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

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

2.     Рассчитывается матрица Якоби значений частных производных в точке начального приближения

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

4.     к вектору начального приближения прибавляется вектор приращений

5.     проверяется условие сходимости и, если оно не достигнуто, то процедура повторяется с п. 2.

 

Метод легко обобщается на систему уравнений любой размерности.   

Для функции F1  n  переменных линейная часть разложения в ряд Тейлора в окрестности точки   записывается так 

После разложения всех уравнений системы и используя введенные ранее обозначения, после преобразования получим систему линейных уравнений порядка  n относительно приращения переменных Δxi

Или, в матричной форме,

В сокращенном виде можно записать так - (F')(Δx) = -(F)  , где матрица значений частных производных – (F') – называется матрицей Якоби или якобианом системы уравнений.

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

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

   ,

где эпсилон – достаточно малое число.

Методы контроля сходимости итерационных методов решения систем

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

1.  Норма (эвклидова или -максимум) вектора невязок  

2.  Эвклидова норма вектора относительных отклонений переменных

  

3.  Норма-максимум вектора относительных отклонений    

Пример

Применим метод Ньютона для решения системы уравнений

Матрица частных производных (в аналитическом виде)

Система линейных уравнений

Может быть решена аналитически или методом Крамера или методом обращения матрицы. Возьмем начальное приближение x = 0,15, y = 0,17

Первая итерация:

Матрица Якоби -

вектор значений функции

Рассчитанный вектор поправок

Новое приближение

 x = 0,15 - 0,0138 = 0,0136196,    y = 0,17 + 0,071604 = 0,241604

Вторая итерация:

Рассчитанный вектор поправок

Новое приближение x=0,168057,    y=0,268688

Третья итерация:

Рассчитанный вектор поправок

Новое приближение x=0,181878,    y=0,282509

На 25-й итерации эвклидова норма вектора невязок составляет  2.77∙10-7, эвклидова норма вектора поправок составляет  1.96∙10-7 и решение сходится к  x = 0.2, y = 0.3 с абсолютной погрешностью менее 1∙10-6 Метод простой итерации при этих же начальных условиях сходится с такой точностью на 33-м шаге, модификация Зейделя – на 31-м шаге.

ТЕМА. Численное решение дифференциальных уравнений

Постановка задачи. Типы задач для ОДУ.

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

u^(p)(x) = f(x, u, u^/, u^//, … u^(p+1)) заменой u^(k)(x) = y_k(x) можно свести к эквивалентной системе из p уравнений первого порядка, представленных в каноническом виде :

y^/_k(x) = y_k+1(x) для 0 ≤ k ≤ p–2 (1)

y^/_p-1(x) = f(x, y₀, y₁, … y_p–1), при этом y₀(x) ≡ u(x). (2)

Покажем такое преобразование на примере уравнения Бесселя: .

Предполагая тождественную замену y₁(x) ≡ y(x) представим систему ОДУ в следующем виде:

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

Известно, что система p-го порядка имеет множество решений, которые в общем случае зависят от p параметров {C₁, C₂, … C_p}. Для определения значений этих параметров, т.е. для выделения единственного решения необходимо наложить p дополнительных условий на u_k(x). По способу задания условий различают три вида задач, для которых доказано существование и единственность решения. Это

1. Задача Коши. Задается координата u_k(x₀) = u_k0 начальной точки интегральной кривой в (p+1) – мерном пространстве (k = 1…p). Решение при этом требуется найти на некотором отрезке x₀ ≤ x ≤ x_max .

2. Краевая задача. Это задача отыскания частного решения системы ОДУ на отрезке a ≤ x ≤ b, в которой дополнительные условия налагаются на значения функции u_k(x) более чем в одной точке этого отрезка.

3. Задача на собственные значения. Кроме искомых функций и их производных в уравнение входят дополнительно m неизвестных параметров λ₁, λ₂, … λ_m, которые называются собственными значениями. Для единственности решения на некотором интервале необходимо задать p+m граничных условий.

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

Метод Эйлера

Рассмотрим задачу Коши для уравнения первого порядка:

u^/(x) = f(x, u(x)), x₀ ≤ x ≤ x_max, u(x₀) = u₀. (3)

В окрестности точки x₀ функцию u(x) разложим в ряд Тейлора:

(4)

И дея этого и последующих методов основывается на том, что если f(x,u) имеет q непрерывных производных, то в разложении можно удержать члены вплоть до O(h^q+1), при этом стоящие в правой части производные можно найти, дифференцируя (3) требуемое число раз. В случае метода Эйлера ограничимся только двумя членами разложения.

Пусть h – малое приращение аргумента. Тогда (4) превратится в

.

Так как в соответствии с (3) u^/(x₀) = f(x₀, u₀), то .

Теперь приближенное решение в точке x₁ = x₀ + h можно вновь рассматривать как начальное условие, т.е. организуется расчет по следующей рекуррентной формуле:

(5),

где y₀ = u₀, а все y_k ­– приближенные значения искомой функции (см. рисунок). В методе Эйлера происходит движение не по интегральной кривой, а по касательной к ней. На каждом шаге касательная находится уже для новой интегральной кривой (что и дало название методу – метод ломаных), таким образом ошибка будет возрастать с отдалением x от x₀.

При приближенное решение сходится к точному равномерно с первым порядком точности. То есть, метод дает весьма низкую точность вычислений: погрешность на элементарном шаге h составляет ½ h² y^//( ½(x_k + x_k+1)), а для всей интегральной кривой порядка h¹. При h = const для оценки апостериорной погрешности может быть применена первая формула Рунге, хотя для работы метода обеспечивать равномерность шага в принципе не требуется.

Метод Эйлера легко обобщается для систем ОДУ. При этом общая схема процесса (5) может быть записана так:

(6),

где i = 1…m – число уравнений, k – номер предыдущей вычисленной точки.

Алгоритм

5. Построить сетку с шагом . Значение шага выбирается из соображений требуемой точности, учитывая, что порядок точности метода Эйлера . Левая граница отрезка, на котором строится сетка, задана начальными условиями задачи Коши.

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

,

где – правая часть дифференциального уравнения.

14. Для проверки сравнить значения приближенного решения со значениями точного решения в узлах сетки.

Модифицированный метод Эйлера

6. Построить сетку с шагом . Значение шага выбирается из соображений требуемой точности, учитывая, что порядок точности модифицированного метода Эйлера . Левая граница отрезка, на котором строится сетка, задана начальными условиями задачи Коши.

15. Решение дифференциального уравнения ищется в виде сеточной функции. Значение известно из начальных условий, все следующие значения рассчитываются по модифицированной формуле Эйлера (формула Эйлера с пересчетом):

,

где – правая часть дифференциального уравнения;

– вычисляется предварительно по формуле простого метода Эйлера.

16. Для проверки сравнить значения приближенного решения со значениями точного решения в узлах сетки.

Метод Рунге–Кутта четвертого порядка точности

Метод Рунге-Кутты II порядка.

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

. (8)

Переходя к приближенному решению y ≈ u и заменяя производные в (8) конечными разностями, получаем в итоге следующее выражение:

, (9)

где 0 ≤ α ≤ 1– свободный параметр. Можно показать, что если f(x,u) непрерывна и ограничена вместе со своими вторыми производными, то решение, полученное по данной схеме, равномерно сходится к точному решению с погрешностью порядка h².

Для параметра α наиболее часто используют следующие значения:

1. α = 1. В этом случае

. Графически это уточнение можно интерпретировать так: сначала по схеме ломаных делается шаг h/2, и находится значение . В найденной точке определяется наклон касательной к интегральной кривой, который и будет определять приращение функции для целого шага, т.е. отрезок [AB] (см. рисунок) будет параллелен касательной, проведенной в точке (x_k + h/2, y(x_k + h/2) ) к интегральной кривой.

2. α = ½. В этом случае

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

Метод Рунге-Кутты IV порядка.

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

, где i = 1…p, p – число уравнений в системе.

; k – номер точки, для которой осуществляется расчет;

;

;

.

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

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

Метод Рунге – Кутта - Мерсона.

Этот метод отличается от метода Рунге – Кутта четвертого порядка возможностью оценивать погрешность на каждом шаге и в зависимости от этого принимать решение об изменении шага. Один из вариантов формул:

;

R_n+1 = 0.2k₄ – 0.3k₃ – 0.1k₅ - погрешность на каждом шаге.

Пусть задана максимальна погрешность . Если , h = h/2 , и n+1 цикл расчета повторяется (с точки x_n, y_n) c новым шагом. Если же

, h = 2h

Автоматический выбор шага позволяет значительно сократить время решения ОДУ.

Схема РКМ обобщается на системы ОДУ аналогично классической схеме Рунге – Кутта.

Метод Адамса.

Метод основан на аппроксимации интерполяционными полиномами правых частей ОДУ.

Пусть с помощью любого из методов, рассмотренных выше, вычислено решение заданного дифференциального уравнения в точках x₁, x₂, x₃ (а в точке x₀ решение и так известно – поставлена задача Коши). Полученные значения функции обозначим как y₀, y₁, y₂, y₃, а значения правой части дифференциального уравнения как f₀, f₁, f₂, f₃, где f_k = f(x_k, y_k). Начиная с четвертой точки, на каждом шаге интегрирования дифференциального уравнения вычисления осуществляются по схеме

P(EC)^{m}E

где P – прогноз решения; Е ­– вычисление f(x,y); С – коррекция решения; m ­– количество итераций коррекции. Схемы такого типа называют «прогноз-коррекция»: это подразумевает сначала приблизительное вычисление решение по формуле низкого порядка, а затем уточнение с учетом полученной информации о поведении интегральной кривой.

Прогноз осуществляется по экстраполяционной формуле Адамса:

. (10)

Коррекция осуществляется по интерполяционной формуле Адамса:

. (11)

Вычисление осуществляется по формуле:

Количество итераций m ≤ p, где p ­– порядок используемого метода. В ходе каждой итерации решается нелинейное уравнение (11) относительно неизвестной y₄ (обычно методом простых итераций).

Иногда в методе Адамса используется схеме PECE на каждом шаге процесса интегрирования, т.е. осуществляется только одна коррекция. В силу сложности вычислений метод используется только в мощных программных пакетах численного анализа. Формулы метода также легко переносятся на решение систем ОДУ первого порядка.

Постановка краевой задачи. Метод стрельбы.

Краевая задача – задача отыскания частного решения системы

(12), 1 ≤ k ≤ p, на отрезке a ≤ x ≤ b, на котором дополнительные условия налагаются на значения функции u_k(x) более чем в одной точке этого отрезка. В качестве примера моно привести задачу нахождения статического прогиба u(x) нагруженной струны с закрепленными концами:

, a ≤ x ≤ b, u(a) = u(b) = 0. Здесь f(x) имеет смысл внешней изгибающей нагрузки на единицу длины струны, деленной на упругость струны.

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

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

(13)

с краевыми условиями

(14)

Сначала выберем некоторое значение , так что . Это уравнение мы можем решить и определить . Таким образом у нас появились два числа и β, которые будут определять задачу Коши , для системы (13), удовлетворять левому краевому условию, но не удовлетворять второму уравнению (14). Тем не менее, решим систему ОДУ с задачей Коши каким-либо из известных нам численных методов. Получим решение вида

(15), зависящее от как от параметра.

Теперь мы должны каким-либо способом менять параметр до тех пор, пока не подберется такое значение, для которого будет выполнено условие (16), т.е. правое краевое условие. Для этого случайным образом берут значения до тех пор, пока среди величин не окажется разных по знаку.

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

Решение уравнений в частных производных.

К уравнениям в частных производных приводят задачи газодинамики, теплопроводности, переноса излучения, электромагнитных полей, процессов переноса в газах, и др. Независимыми переменными в физических задачах обычно являются время t, координаты , скорости частиц . Пример – уравнение теплопроводности

, (17)

где U ­– температура, – теплоемкость, – коэффициент теплопроводности, q – плотность источников тепла.

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

Рассмотрим в качестве примера одномерную задачу, близкую по смыслу к (17):

. (18)

Здесь 0 ≤ x ≤ a, 0 ≤ t ≤ T.

Граничные условия:

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

Введем равномерную прямоугольную сетку по x и t с шагом h и τ соответственно и заменить производные соответствующими разностными отношениями. Тогда

. (19)

Здесь 1 ≤ k ≤ N-1 – число точек по координате x; 0 ≤ m ≤ M – число точек по координате t. Число неизвестных в (19) больше числа уравнений, недостающие уравнения выводятся из начальных и граничных условий:

; 0 ≤ k ≤ N.

; ; 0 ≤ m ≤ M.

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

, где 1 ≤ k ≤ N-1.

; . (20)

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

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

Алгоритм метода РК 4 порядка

7. Построить сетку с шагом . Значение шага выбирается из соображений требуемой точности, учитывая, что порядок точности метода . Левая граница отрезка, на котором строится сетка, задана начальными условиями задачи Коши.

17. Решение дифференциального уравнения ищется в виде сеточной функции. Значение известно из начальных условий, все следующие значения рассчитываются по формулам Рунге–Кутта:

,

,

,

,

,

где – правая часть дифференциального уравнения;

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

18. Для проверки сравнить значения приближенного решения со значениями точного решения в узлах сетки.