Тема 7. Решение систем нелинейных алгебраических уравнений

Определения

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

( 7.1)

где f_i (i=1, 2, …, n ) - функции вещественных переменных x₁, x₂, …, x_n. Обозначим упорядоченную совокупность n чисел x = ( x₁, x₂, …, x_n)^TÎ H и F(x) = (f(x₁), f( x₂), …, f(x_n))^T. Тогда система уравнений(7.1) в некотором линейном пространстве H размерности n запишется в операторном виде: F(x)=0, где F: H® H нелинейное отображение. Такие системы решают итерационными методами.

Многие одношаговые итерационные методы для решения системы F(x)=0 можно записать в канонической форме:

. (7.2)

где k - номер итерации x^k = ( x^k₁, x^k₂, …, x^k_n)^T, t_k₊₁ - числовой параметр, - матрица n x n, имеющая обратную матрицу.

Метод (7.2) называется явным, если B_k₊₁ =E для всех k, и неявным в противном случае.

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

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

Преобразуем эту систему к виду:

Построим итерационный процесс

(7.3)

Итерационный метод (7.3) сходится, если

a) в области R={a≤ x ≤ b, c≤ y ≤ d} существует единственное решение;

б) в области R выполнены условия

или

Начальное приближение (x⁰, y⁰) должно принадлежать области R.

Метод релаксации

Если B_k₊₁ =E, t_k₊₁=t, это стационарный итерационный метод, который можно записать в виде: x^k=S(x^k), где S(x)=x-t F(x). Метод сходится, если . В данном случае , а F΄ (x) - матрица Якоби.

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

Пусть известно приближение на k -м шаге x^k= . Выпишем разложение функции f_i(x₁, x₂, …, x_n) по формуле Тейлора в точке x^k до второго порядка:

Тогда система (7.1) заменится системой уравнений:

i=1, 2, …, n, ( 7.4)

линейной относительно приращений x_j-x_j^k, j=1, 2, …n. Решение x=(x_1,x₂, …, x_n)^T системы (7.3) примем за следующее итерационное приближение и обозначим x^k⁺¹ = .

Таким образом, итерационный метод Ньютона для системы (7.1) определяется системой уравнений:

i=1, 2, …, n, (7.5)

из которой последовательно, начиная с заданного x⁰=( x₁⁰, x₂⁰, …, x_n⁰)^T, находятся векторы x^k, k=1, 2, ….. Систему (7.4) можно записать в векторном виде:

k=1, 2, …,

где x⁰- заданный вектор, F΄ (x) - матрица Якоби.

(7.6)

Если обозначить z^k=x^k^{+ 1}- x^k и если (F΄ (x^k))^-1 существует, то, решив систему линейных алгебраических уравнений:

F΄ (x^k) z^k= - F(x^k),

(k+1) -е приближение найдем из формулы

x^k+1=z^k + x^k.

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

или .

Задания

1. Найти решение системы

- методом Ньютона

- методом простой итерации с заданной точностью ε . Начальное приближение (x⁰, y⁰) найти графически.

Варианты

1.	7.
2.	8.
3.	9.
4.	10.
5.	11.
6.	12.

Тема 8. Решение обыкновенных дифференциальных уравнений

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

=f(x, y), y(a)=y₀; xÎ [a, b]. (8.1)

Метод численного интегрирования

Интегрируя (8.1), заменим эту задачу эквивалентным ей интегральным уравнением типа Вольтера:

y(x)= y₀+ f(t, y(t)) dt. (8.2)

Решим уравнение (8.2) методом последовательных приближений. В результате получим итерационный процесс

y^(k)(x) = y₀+ f(t, y^(k-1)(t)) dt, (8.3)

Зададим сетку x_i=a+ih; i=0, …, m, где h= - шаг сетки. Из уравнения (8.3) имеем

y^(k)(x_i) = y_i^(k) = y₀+ f(t, y^(k-1)(t)) dt, (8.4)

где y⁽⁰⁾(t)=y₀.

Интеграл в правой части (8.4) вычисляется с помощью численной квадратуры

y_i^(k) –y₀ = h A_j⁽ⁱ⁾ f(x_j, y_j^(k-1)). (8.5)

Коэффициенты A_j⁽ⁱ⁾ системы (8.5) находятся из решения системы линейных алгебраических уравнений:

= A_j⁽ⁱ⁾ j^k^-1; k=1, …n+1; i=1, …n. (8.6)

Система (8.6) получена из условия, что формула (8.5) точна для всех функций вида x, x², x³, …, xⁿ⁺¹в точках x_i = ih, i=0, …, n.

Пусть n=2, дано три точки (x₀, x₁, x₂). Тогда матрица коэффициентов A_j⁽ⁱ⁾ равна

A= .

Для m=4, пять точек (x₀, x₁, x₂, x₃, x₄).

A= .

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

Шаг 1. Исходный интервал делим на 4 части, т. е. задаем m = 4. Определяем сетку x_i=a+ih₀, i=0, 1, …, 4, где h₀= (a+b)/4 - шаг исходной сетки. Вычисляем по формулам (8.5) y_i⁽ⁿ⁾, i=1, …, 4, пока не выполнится условие max < e, по всем i, где e - заданная точность.

Шаг 2.Делим интервал на 8 частей, т. е. уменьшаем шаг вдвое и находим y_i на сетках:

=a+ih₁, i=1, …, 4; h₁=h₀/2;

=(b+a)/2+ih₁, i=1, …, 4.

Затем проверяем условие max < e в узлах первой сетки (с шагом h₀). Если оно выполняется, то вычисления заканчиваем. Если не выполняется, то исходный интервал уменьшаем вдвое, т. е. a₁=a; b₁=(a+b)/2, и переходим на шаг 1, где a=a₁; b=b₁.

Метод Эйлера

Для решения задачи Коши (8.1) составляют таблицу значений y_k=y(x_k), где x_k=x₀+kh ( k=0, 1, …, n), h=(b-a)/n, y(x₀)=y₀, x₀=a, [a, b] – отрезок, на котором ищется решение. Значение y_k₊₁ определяется по формуле

y_k₊₁=y_k + h f(x_k, y_k), k=0, 1, …, n -1, y(x₀)=y₀. (8.7)

Погрешность вычислений на каждом шаге составляет

R_k=0.5h² (e), где x_k£ e £ x_k₊₁.

Модифицированный метод Эйлера с уточнением состоит в следующем: сначала вычисляют

y⁽⁰⁾_k+1= y_k+ h f(x_k, y_k), k=0, 1, …n -1,

а затем это значение уточняют по формуле

y⁽ⁱ⁾_k₊₁=y_k+ [f(x_k, y_k)+f(x_k₊₁, y⁽ⁱ^-1)_k₊₁)], где i=0, 1, 2, … - номер итерации. (8.8)

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

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

Для решения задачи Коши (8.1) семейство методов Рунге-Кутта описывается следующим выражением:

y(x+h)=y(x)+ , (8.9)

где h – шаг сетки; число q называется порядком точности метода (8.9).

k₁=h f(x, y),

k₂=h f(x+a₂h, y+b₂₁k₁),

…

k_i=h f(x+a_Ih, y+b_i1k₁+b_i2k₂+…+b_ii-1k_i-1); i=1, …, q.

Здесь p_i, a_i, b_ik – коэффициенты.

При q = 1 имеем метод первого порядка точности:

y_i+1=y_i+k₁,

k₁=h f(x_i, y_i),

y₀=y(x₀).

Метод второго порядка точности (при q = 2) имеет вид:

y_i+1=y_i+ (k₁+ k₂),

k₁=h f(x_i, y_i),

k₂=h f(x_i+h, y_i+k₁),

y₀=y(x₀).

Наиболее распространен на практике метод четвертого порядка точности (при q=4):

y_i+1=y_i+ (k₁+2(k₂+k₃)+k₄),

где k₁=h f(x_i, y_i),

k₂= h f(x_i+h/2, y_i+ k₁ /2),

k₃= h f(x_i+h/2, y_i+ k₂ /2),

k₄= h f(x_i+ h, y_i+ k₃) .

Для нахождения решения с заданной точностью e, численная реализация методов Рунге-Кутта выполняется следующим образом: задается сетка x_i=a+ih₀, i =1, …, n; h₀=(b-a)/n, далее на каждом I -м шаге в точке x=x_i вычисляют два значения y_i⁽¹⁾ и y_i⁽²⁾:

y_i⁽¹⁾=y_i-1+ p_j(h₀)k_j(h₀, x_i-1, y_i-1),

=y_i-1+ p_j(h₁)k_j(h₁, x_i-1, y_i-1),

y_i⁽²⁾= + p_j(h₁)k_j(h₁, x_i-1, ), где h₁=h₀/2.

Если выполнено условие | y_i⁽¹⁾- y_i⁽²⁾| < e, где e - заданная точность, то следующее y_i₊₁ вычисляется c тем же шагом h₀. В противном случае полагают h₀=h₁, h₁=h₀ /2.

Метод Адамса

В вышеизложенных методах для задачи Коши значение y_i₊₁ зависело только от информации в предыдущей точке x_i (одношаговые методы). В многошаговых методах используют информацию в нескольких предыдущих точках x_i, x_i_-1, x_i_-2, …..

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

y_i+1⁽^р⁾=y_i+ (55 f_i- 59 f_i-1+ 37 f_i_-2- 9 f_i-3),

f_i+1^(p)=f(x_i+1, y_i+1^(p) ),

y_i+1=y_i+ (9 f_i+1^(p)+19 f_i - 5 f_i-1+ f_i_-2), i=3, …n.

Вычисленное значение y_i₊₁^(р), являющееся «прогнозом» для y_i₊₁, затем y_i₊₁^(р)используют для вычисления приближенного значения f_i₊₁, которое, в свою очередь, используют в формуле для вычисления y_i₊₁. В начале работы необходимо вычислить значения в точках y_i (i=1, 2, 3) с помощью одношагового метода того же порядка точности.

Метод Милна

Пусть для решения уравнения =f(x, y), кроме начального условия y(x₀)=y₀, найдены значения искомой функции y(x_i)=y_i в точках x_i=x₀+ih (i=1, 2, 3). Последующие значения y_i при i=4, 5, … определим, используя формулы метода Милна:

y_i^пред=y_i_-4+ (2 f_i_-3- f_i_-2+ 2 f_i_-1) (прогноз),

y_i^кор=y_i_-2+ (f_i_-2– 4 f_i_-1+f_i^пред), где f_i^пред=f(x_i, y_i^пред ) (коррекция).

Абсолютная погрешность e_i значения y_i^кор приближенно определяется по формуле e_i@ | y_i^кор- y_i^пред|. Если точность результата достаточна, то полагают y_i @ y_i^кор.