Реализация некоторых численных методов

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

В ряде случаев весьма удобным приёмом уточнения корня уравнения является метод последовательных приближений (метод итераций).

Пусть с точностью $$\varepsilon$$ необходимо найти корень уравнения $$f(x) = 0$$, принадлежащий интервалу изоляции $$[a,b]$$. Функция $$f(x)$$ и ее первая производная непрерывны на этом отрезке.

Для применения этого метода исходное уравнение $$f(x) = 0$$ должно быть приведено к виду $$x=\varphi(x)$$.

В качестве начального приближения может быть выбрана любая точка интервала $$[a,b]$$.

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

В результате итерационный процесс поиска реализуется рекуррентной формулой. Процесс поиска прекращается, как только выполняется условие $$\lvert x_n-x_{n-1} \rvert \le \varepsilon$$ или число итераций превысит заданное число $$N$$.

Для того, чтобы последовательность $$x_1,x_2,... ,x_n$$ приближалась к искомому корню, необходимо, чтобы выполнялось условие сходимости $$\lvert \varphi'(x) \rvert <1$$.

Пример реализации метода итераций представлен ниже:

(%i1)	f:exp(-x)-x$
beta:0.1$ x1:1$ x0:0$ eps:0.000001$ p:0$
while abs(x1-x0)>eps do
		(x0:x1, p:p+1, x1:float(x0+beta*(subst(x0,x,f))))$
print("Число итераций ",p," ","Решение ",float(x1),
" Невязка ",float(abs(x1-x0)))$

8.1.3 Метод Ньютона (метод касательных)

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

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

Рассмотрим нелинейное уравнение $$f(x) = 0$$, для которого необходимо найти корень на интервале $$[a,b]$$ с точностью $$\varepsilon$$.

Метод Ньютона основан на замене исходной функции $$f(x)$$, на каждом шаге поиска касательной, проведённой к этой функции. Пересечение касательной с осью $$X$$ дает приближение корня.

Выберем начальную точку $$x_0 = b$$ (конец интервала изоляции). Находим значение функции в этой точке и проводим к ней касательную, пересечение которой с осью $$X$$ дает первое приближение корня $$x_1$$:

Поэтому $$x_1=x_0-\dfrac{f(x_0)}{f'(x_0)}$$.

В результате, итерационный процесс схождения к корню реализуется рекуррентной формулой

Процесс поиска продолжаем до тех пор, пока не выполнится условие: $$\lvert x_{n+1}-x_n \rvert \le \varepsilon$$, откуда $$\lvert \cfrac{f(x_n)}{f'(x_n)}\rvert \le\varepsilon$$.

Метод обеспечивает быструю сходимость, если выполняется условие: $$f(x_0)\cdot f''(x_0) > 0$$, т.е. первую касательную рекомендуется проводить в той точке интервала $$[a,b]$$, где знаки функции $$f(x_0)$$ и ее кривизны $$f''(x_0)$$ совпадают.

Пример реализации метода Ньютона в Maxima представлен ниже:

(%i1)	newton(f,x0,eps):=block([df,xn,xn0,r,p],
	xn0:x0, df:diff(f,x),
	p:0, r:1,
	while abs(r)>eps do (
		p:p+1, xn:xn0-float(subst(xn0,x,f)/subst(xn0,x,df)),
		print("x0,x1 ",xn0,xn),r:xn-xn0, xn0:xn
	),
	[xn,p])$

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

(%i2)	f:exp(-x)-x$
	eps:0.000001$ xrez:newton(f,1,eps)$
	print("Решение ",xrez[1]," Число итераций ",xrez[2],
	" Невязка ",subst(xrez[1],x,f))$

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

8.1.4 Модифицированный метод Ньютона (метод секущих)

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

Тогда рекуррентная формула метода Ньютона приобретёт вид: $$x_{n+1}=x_n-\frac{f(x_n)}{f'(x_n)}=x_n-\frac {f(x_n) \Delta x}{\Delta f(x_n)}=\\ =x_n-\frac {f(x_n) \Delta x}{f(x_n+ \Delta x)-f(x_n)},$$ где $$\Delta x \approx \varepsilon$$.

Реализация данного метода в Maxima представлено ниже:

(%i1)	secant(f,sp,eps):=block([x0,x1,d,y,r],
	x0:sp[1],x1:sp[2],
	p:0, r:x1-x0, d:float(subst(x0,x,f)),
	while abs(r)>eps do (
		p:p+1, y:float(subst(x1,x,f)), r:r/(d-y)*y,
		d:y, x1:x1+r
	),
	x1)$
(%i2)	f:exp(-x)-x$
	eps:0.000001$ xrez:secant(f,[-2,2],eps)$
	print("Решение ",xrez," Невязка ",subst(xrez,x,f))$

Особенности программирования для Maxima, использованные в этом примере, аналогичны приведённым выше в примере, касающемся метода половинного деления.

8.1.5 Метод хорд

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

При этом в процессе поиска семейство хорд может строится:

В результате итерационный процесс схождения к корню реализуется рекуррентной формулой:

Процесс поиска продолжается до тех пор, пока не выполнится условие

Метод обеспечивает быструю сходимость, если $$f(z) \cdot f"(z) > 0$$, т.е. хорды фиксируются в том конце интервала $$[a,b]$$, где знаки функции $$f(z)$$ и ее кривизны $$f''(z)$$ совпадают.

8.2 Численное интегрирование

Задача численного интегрирования состоит в замене исходной подинтегральной функции $$f(x)$$, для которой трудно или невозможно записать первообразную в аналитике, некоторой аппроксимирующей функцией $$\varphi(x)$$. Такой функцией обычно является полином (кусочный полином) $$\varphi (x) = \sum\limits_{i = 1}^n {{c_i}{\phi _i}(x)}$$.

Таким образом $$\text{I}} = \int\limits_a^b {f(x)dx = \int\limits_a^b {\phi (x)dx + R},$$ где $$R = \int\limits_a^b {r(x)dx}$$ — априорная погрешность метода на интервале интегрирования, а $$r(x)$$ — априорная погрешность метода на отдельном шаге интегрирования.

8.2.1 Обзор методов интегрирования

Методы вычисления однократных интегралов называются квадратурными (для кратных интегралов — кубатурными), и делятся на следующие группы:

8.2.2 Метод прямоугольников

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

Формулы метода прямоугольников можно получить из анализа разложения функции $$f(x)$$ в ряд Тейлора вблизи некоторой точки $$x = x_i$$:

Рассмотрим диапазон интегрирования от $$x_i$$ до $$x_i + h$$, где $$h$$ — шаг интегрирования. Вычислим интеграл от исследуемой функции на этом промежутке: $$\int\limits_{x_i}^{{x_i} + h} f(x)dx = \left. {x \cdot f({x_i})} \right|_{{x_i}}^{{x_i} + h} + \frac{{{{\left( {x - {x_i}} \right)}^2}}}{2}\left. {f'({x_i})} \right|_{{x_i}}^{{x_i} + h} + \\ +\frac{{{{\left( {x - {x_i}} \right)}^3}}}{{3 \cdot 2!}}\left. {f''({x_i})} \right|_{{x_i}}^{{x_i} + h} +\dots = \\ =f({x_i})h + \frac{{{h^2}}}{2}f'({x_i}) + O({h^3})=f({x_i})h + {r_i}.$$ таким образом, на базе анализа ряда Тейлора получена формула правых (или левых) прямоугольников и априорная оценка погрешности $$r$$ на отдельном шаге интегрирования. Основной критерий, по которому судят о точности алгоритма — степень при величине шага в формуле априорной оценки погрешности. В случае равного шага $$h$$ на всем диапазоне интегрирования общая формула имеет вид $$\int\limits_a^b {f(x)dx = h\sum\limits_{i = 0}^{n - 1} {f({x_i}) + R} },$$ где n — число разбиений интервала интегрирования, $$R = \sum\limits_{i = 0}^{n - 1} {{r_i}} = \cfrac{h}{2} \cdot h\sum\limits_{i = 0}^{n - 1} {f'({x_i})} = \cfrac{h}{2}\int\limits_a^b {f'(x)dx}$$. Полученная оценка справедлива при наличии непрерывной производной подинтегральной функции $$f'(x)$$.

8.2.3 Метод средних прямоугольников

Здесь на каждом интервале значение функции считается в средней точке отрезка $$[x_i,x_i + h]$$, то есть $$\int\limits_{{x_i}}^{{x_i} + h} {f(x)dx =hf(\bar x) + {r_i}}$$.

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

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

(%i1)	intpr(f,n,a,b):=block([h,i,s,_x],h:(b-a)/n, _x:a+h/2, s:0,
	for i:1 thru n do (s:s+float(subst(_x,x,f)),_x:_x+h),s:s*h)$
(%i2)	intpr(x^2,100,-1,1);

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

Аппроксимация в этом методе осуществляется полиномом первой степени. На единичном интервале

В случае равномерной сетки ($$h$$ = const)

При этом $${r_i} = - \cfrac{{{h^3}}}{{12}}f''({x_i})$$, а $$R = -\cfrac{{{h^3}}}{{12}}\int\limits_a^b {f''(x)} dx$$.

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

Пример реализации метода трапеций в виде функции приведен ниже:

(%i1)	f:x^2;

(%i2)	inttrap(f,n,a,b):=block([h,i,s],h:(b-a)/n,
	s:(float(subst(a,x,f))+float(subst(b,x,f)))/2,
	for i:1 thru n-1 do (s:s+float(subst(a+i*h,x,f))), s:s*h)$
(%i3)	inttrap(f,100,-1,1);

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

(%i4)	inttrap(x*sin(x),100,-1,1);

8.2.5 Метод Симпсона

При использовании данного метода подинтегральная функция $$f(x)$$ заменяется интерполяционным полиномом второй степени $$P(x)$$ — параболой, проходящей через три соседних узла. Рассмотрим два шага интегрирования ($$h = const = x_{i+1} - x_i$$), то есть три узла $$x_0, x_1,x_2$$, через которые проведем параболу, воспользовавшись уравнением Ньютона:

Пусть $$z = x - x_0$$, тогда

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

В итоге $$\int\limits_{{x_0}}^{{x_2}} {f(x)dx} = \cfrac{h}{3} \left({f_0} +4{f_1} + {f_2} \right) + r$$.

Для равномерной сетки и чётного числа шагов n формула Симпсона принимает вид: $$\int\limits_a^b {f(x)dx} = \frac{h}{3} \left( {f_0} + 4{f_1} + 2{f_2} + 4{f_3} + \dots + 2{f_{n - 2}} + 4{f_{n - 1}} + {f_n}\right) + R,$$ где $$r = - \cfrac{h^5}{90}{f^{IV}}({x_i})$$, а $$R = -\cfrac{h^4}{180}\int\limits_a^b {f^{IV}(x)dx}$$ в предположении непрерывности четвёртой производной подинтегральной функции.

Реализация метода Симпсона средствами Maxima представлена в следующем примере:

(%i1)	intsimp(f,n,a,b):=block([h,h2,i,_x,s],h:(b-a)/n, h2:h/2,
	s:(float(subst(a,x,f))+float(subst(b,x,f)))/2+
				2*float(subst(a+h2,x,f)),
	_x:a,
	for i:1 thru n-1 do (_x:_x+h,
		s:s+2*float(subst(_x+h2,x,f))+float(subst(_x,x,f))),
	s:s*h/3)$
(%i2)	intsimp(x^2,100,-1,1);

8.3 Методы решения систем линейных уравнений

8.3.1 Общая характеристика и классификация методов решения

Рассмотрим систему линейных алгебраических уравнений (сокращенно — СЛАУ): $$A \cdot \bar x = \bar f,$$ где $$A$$ — матрица $$m \times m,\ \bar x = {({x_1},{x_2},\dots,{x_m})^T}$$ — искомый вектор, $$\bar f = {({f_1},{f_2},\dots,{f_m})^T}$$ — заданный вектор. Будем предполагать, что определитель матрицы A отличен от нуля, т.е. решение системы (8.1) существует.

Методы численного решения системы (8.1) делятся на две группы: прямые методы ("точные") и итерационные методы.

Прямыми методами называются методы, позволяющие получить решение системы (8.1) за конечное число арифметических операций. К этим методам относятся метод Крамера, метод Гаусса, LU-метод и т.д.

Итерационные методы (методы последовательных приближений) состоят в том, что решение системы (8.1) находится как предел последовательных приближений $${\bar x^{(n)}}$$ при $$n \to \infty$$, где $$n$$ — номер итерации. При использовании методов итерации обычно задается некоторое малое число $$\varepsilon$$ и вычисления проводятся до тех пор, пока не будет выполнена оценка $$\left\| {{{\bar x}^{(n)}} - \bar x} \right\| < \varepsilon$$. К этим методам относятся метод Зейделя, Якоби, метод верхних релаксаций и т.д.

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

8.3.2 Метод Гаусса

Запишем систему (8.1), в развернутом виде

Метод Гаусса состоит в последовательном исключении неизвестных из этой системы. Предположим, что $${a_{11}} \ne 0$$. Последовательно умножая первое уравнение на $$- \cfrac{{{a_{i1}}}}{{{a_{11}}}}$$ и складывая с $$i$$-м уравнением, исключим $$x_1$$ из всех уравнений кроме первого. Получим систему $$\begin{aligned} a_{11}{x_1} + a_{12}{x_2} + \dots + a_{1m}{x_m} = {f_1}, \\ a_{22}^{(1)}{x_2} + \dots + a_{2m}^{(1)}{x_m} = f_2^{(1)}, \\ \dots\dots\dots\dots\dots\dots\dots\dots\dots\dots \dots\dots\dots \\ a_{m2}^{(1)}{x_2} + \dots + a_{mm}^{(1)}{x_m} = f_m^{(1)}. \end{aligned}$$ где $$a_{ij}^{(1)} = {a_{ij}} - \frac{a_{i1}a_{1j}}{a_{11}}, \quad f_i^{(1)} = f_i - \frac{a_{i1}f_1}{a_{11}}, \quad i,j = 2,3, \dots, m.$$

Аналогичным образом из полученной системы исключим $$x_2$$. Последовательно, исключая все неизвестные, получим систему треугольного вида $$\begin{aligned} {a_{11}}{x_1} + {a_{12}}{x_2} + \dots + {a_{1k}}{x_k} + \dots +{a_{1m}}{x_m} = f_1, \\ a_{22}^{(1)}{x_2} + \dots + a_{2k}^{(1)}{x_k} + \dots + a_{2m}^{(1)}{x_m} = f_2^{(1)}, \\ \dots\dots\dots\dots\dots\dots\dots\dots\dots\dots\dots\dots \dots\dots\dots, \\ a_{m - 1,m - 1}^{(m - 1)}{x_{m - 1}} + a_{m - 1,m}^{(m - 1)}{x_m} = f_{m - 1}^{(m - 1)},\\ a_{m,m}^{(m - 1)}{x_m} = f_m^{(m - 1)}.\\ \end{aligned}$$

Описанная процедура называется прямым ходом метода Гаусса. Заметим, что ее выполнение было возможно при условии, что все $$a_{i,i}^{(l)},\ l = 1,2,\dots ,m-1$$ не равны нулю. Выполняя последовательные подстановки в последней системе, (начиная с последнего уравнения) можно получить все значения неизвестных.

Эта процедура получила название обратный ход метода Гаусса.

Метод Гаусса может быть легко реализован на компьютере. При выполнении вычислений, как правило, промежуточные значения матрицы A не представляют интереса. Поэтому численная реализация метода сводится к преобразованию элементов массива размерности $$m \times (m+1)$$, где $$m + 1$$-й столбец содержит элементы правой части системы.

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

Один из основных недостатков метода Гаусса связан с тем, что при его реализации накапливается вычислительная погрешность.

Для систем порядка m число действий умножения и деления близко $$\cfrac{m^3}{3}$$ и быстро растет с величиной $$m$$.

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

Пример реализации метода Гаусса в Maxima приведен в функции ниже (применен метод без выбора главного элемента):

(%i1)	gauss(a0,b0,n):=block([a,b,i,j,k,d],
	a:copymatrix (a0), b:copymatrix(b0), x:copymatrix(b0),
	for i:1 thru n-1 do
	(
		for k:i+1 thru n do
		(
			d:a[k,i]/a[i,i],
			for j:i+1 thru n do (a[k,j]:a[k,j]-a[i,j]*d),
			b[k,1]:b[k,1]-b[i,1]*d
		)
	),
	for i:n thru 1 step -1 do
	(
		for j:i+1 thru n do (
			b[i,1]:b[i,1]-a[i,j]*x[j,1]),
		x[i,1]:b[i,1]/a[i,i]
	),
	x)$

Пример обращения к функции, реализующей метод Гаусса:

(%i2)	aa:matrix([3,1,1],[1,3,1],[1,1,3]); bb:matrix([6],[6],[8]);
	zz:gauss(aa,bb,3);

Проверка вычислений показывает, что перемножение матрицы $$A$$ на вектор решения $$zz$$ дает вектор, совпадающий с вектором правых частей $$bb$$.

(%i5)	aa.zz;

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

Выполняемые в методе Гаусса преобразования прямого хода, приведшие матрицу $$A$$ системы к треугольному виду позволяют вычислить определитель матрицы

Метод Гаусса позволяет найти и обратную матрицу. Для этого необходимо решить матричное уравнение $$A \cdot X = E,$$ где $$E$$ – единичная матрица. Его решение сводится к решению $$m$$ систем $$A{\bar x^{(j)}} = {\bar \delta ^{(j)}},\begin{array}{*{20}{c}} {}{j = 1,2,\dots,m,}$$ у вектора $${\bar \delta ^{(j)}}\ j$$–я компонента равна единице, а остальные компоненты равны нулю.

8.3.3 Метод квадратного корня

Метод квадратного корня основан на разложении матрицы $$A$$ в произведение $$A=S^{T}S,$$ где $$S$$ — верхняя треугольная матрица с положительными элементами на главной диагонали, $$S^T$$ — транспонированная к ней матрица.

Пусть $$A$$ — матрица размером $$m\times m$$. Тогда

Из условия (8.2) получаются уравнения

Так как матрица $$A$$ симметричная, не ограничивая общности, можно считать, что в системе (8.3) выполняется неравенство $$i\le j$$. Тогда (8.3) можно переписать в виде

В частности, при $$i = j$$ получится

Далее, при $$i < j$$ получим

По приведённым формулам находятся рекуррентно все ненулевые элементы матрицы $$S$$.

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

Решения этих систем находятся по рекуррентным формулам:

Всего метод квадратного корня при факторизации $$A = S^TS$$ требует примерно $$\cfrac{m^3}{3}$$ операций умножения и деления и m операций извлечения квадратного корня.

Пример функции, реализующей метод квадратного корня:

(%i1)	holetsk(a0,b0):=block([L,Lt,x,y,i,j,k,n],
	n:length(a0), L:zeromatrix(n,n),
	for i:1 thru n do (
		for j:1 thru i-1 do (
		s:0,
		for k:1 thru j-1 do (s:s+L[i,k]*L[j,k]),
		L[i,j]:1/L[j,j]*(a0[i,j]-s)
		),
		s:0,
		for k:1 thru i-1 do (s:s+L[i,k]^2),
		L[i,i]:sqrt(a0[i,i]-s)
		),Lt:transpose(L),
	y:zeromatrix(n,1), x:zeromatrix(n,1),
	for i:1 thru n do (
		s:0,
		for k:1 thru i-1 do (s:s+L[i,k]*y[k,1]),
			y[i,1]:(b0[i,1]-s)/L[i,i]
	),
	for i:n thru 1 step -1 do (
		s:0,
		for k:n thru i+1 step -1 do (s:s+Lt[i,k]*x[k,1]),
		x[i,1]:(y[i,1]-s)/Lt[i,i]
	),x
	)$

Тест данной функции (решение системы $$Ax = B, B$$ — матрица $$n\times 1,\ A$$ — квадратная симметричная матрица $$n \times n$$, результат решения — вектор $$n \times 1$$):

(%i2)	A:matrix([4,1,1],[1,4,1],[1,1,4])$ B:matrix([1],[1],[1])$

Результаты вычислений:

(%i4)	x:holetsk(A,B);

Проверка:

(%i5)	A.x;

8.3.4 Корректность постановки задачи и понятие обусловленности

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

Говорят, что задача поставлена корректно, если решение существует и единственно и если оно непрерывно зависит от входных данных. Последнее свойство называется также устойчивостью относительно входных данных.

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

Известно, что решение задачи (8.1) существует тогда и только тогда, когда $$\det A \ne 0$$. В этом случае можно определить обратную матрицу $${A^{ - 1}}$$ и решение записать в виде $$\bar x = {A^{ - 1}}\bar f$$.

Исследование устойчивость задачи (8.1) сводится к исследованию зависимости ее решения от правых частей $$\bar f$$ и элементов $$a_{ij}$$ матрицы $$A$$. Для того чтобы можно было говорить о непрерывной зависимости вектора решений от некоторых параметров, необходимо на множестве $$m$$-мерных векторов принадлежащих линейному пространству $$\mathbb H$$, ввести метрику.

В линейной алгебре предлагается определение множества метрик $${l_p}$$ — норма $$\left\|{\bar x} \right\|_p = \left( \sum\limits_{i = 1}^m \left| x_i \right|^p \right)^{1/p}$$ из которого легко получить наиболее часто используемые метрики

Подчиненные нормы матриц определяемые как $$\left\| A \right\| = \mathop {\sup }\limits_{0 \ne x \in \mathbb H} \cfrac{{\left\| {A\bar x} \right\|}}{{\left\| {\bar x} \right\|}}$$, соответственно записываются в следующем виде:

Обычно рассматривают два вида устойчивости решения системы (8.1):

Наряду с исходной системой (8.1) рассмотрим систему с "возмущенными" правыми частями $$A \cdot \tilde x = \tilde f,$$ где $$\tilde f = \bar f + \delta \bar f$$ возмущенная правая часть системы, а $$\tilde x = \bar x + \delta \bar x$$ возмущенное решение.

Можно получить оценку, выражающую зависимость относительной погрешности решения от относительной погрешности правых частей $$\frac{{\left\| {\delta \bar x} \right\|}}{{\left\| {\bar x} \right\|}} \leqslant {M_A}\frac{{\left\| {\delta \bar f} \right\|}}{{\left\| {\bar f} \right\|}},$$ где $${M_A} = \left\| A \right\| \cdot \left\| {{A^{ - 1}}} \right\|$$ — число обусловленности матрицы $$A$$ (в современной литературе это число обозначают как $$cond(A)$$). Если число обусловленности велико ($${M_A}\sim{10^k}, \,\,\,k > 2$$), то говорят, что матрица $$A$$ плохо обусловлена. В этом случае малые возмущения правых частей системы (8.1), вызванные либо неточностью задания исходных данных, либо вызванные погрешностями вычисления существенно влияют на решение системы.

Если возмущение внесено в матрицу $$A$$, то для относительных возмущений решения имеет место следующая оценка:

В Maxima матричные нормы вычисляются посредством функции $$mat_norm$$. Синтаксис вызова: $$mat_norm(M,type)$$, где $$M$$ — матрица, $$type$$ — тип нормы, $$type$$ может быть равен 1 (норма $${\left\| A \right\|_1}$$), $$inf$$ (норма $${\left\| A \right\|_\infty }$$), $$frobenius$$ (норма $$\left\| A \right\|_2$$).

Пример вычисления указанных видов нормы в Maxima:

(%i1)	A:matrix([1,2,3],[4,5,6],[7,8,9]);

(%i2)	mat_norm(A,1);

(%i3)	mat_norm(A,inf);

(%i4)	mat_norm(A,frobenius);

Вычислим число обусловленности для плохо и хорошо обусловленных матриц:

(%i1)	A:matrix([1,1],[0.99,1]);

(%i2)	nrA:mat_norm(A,frobenius);

(%i3)	A1:invert(A);

(%i4)	nrA1:mat_norm(A1,frobenius);

(%i5)	MA:nrA*nrA1;

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

Аналогичным путём (с использованием нормы Фробениуса) для матрицы $$B=\left[ \begin{matrix} 1 0 \\ 0 1 \end{matrix} \right]$$ число обусловленности составило $$MB = 2$$.

8.3.5 О вычислительных затратах

Один из важных факторов предопределяющих выбор того или иного метода при решении конкретных задач, является вычислительная эффективность метода. Учитывая, что операция сложения выполняется намного быстрее, чем операция умножения и деления, обычно ограничиваются подсчетом последних. Для решения СЛАУ методом Гаусса без выбора главного элемента требуется $$\cfrac{{{m^3}}}{3} + {m^2} - \cfrac{m}{3}$$ умножений и делений, решение СЛАУ методом квадратного корня требует $$\cfrac{{{m^3}}}{6} + \cfrac{3m^2}{2} +\cfrac{m}{3}$$ и $$m$$ операций извлечения корней. При больших значениях размерности $$m$$, можно сказать, что вычислительные затраты на операции умножения и деления в методе Гаусса составляют величину $$O\left( {\cfrac{{{m^3}}}{3}}\right)$$, в методе квадратных корней $$O\left( {\cfrac{{{m^3}}}{6}} \right)$$.

8.4 Итерационные методы

В приближенных или итерационных методах решение системы линейных алгебраических уравнений является пределом итерационной последовательности, получаемой с помощью этих методов. К ним относятся: метод простой итерации, метод Зейделя и др. Итерационные методы выгодны для системы специального вида, со слабо заполненной матрицей очень большого вида порядка $$10^3 \dots 10^5$$.

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

8.4.1 Матричная формулировка итерационных методов решения систем линейных уравнений

При использовании СКМ Maxima вполне обосновано использование и матричной формулировки итерационных методов.

Рассмотрим решение системы $$Ax = f$$ ($$A$$ — квадратная матрица, $$f$$ — вектор правых частей, $$x$$ — вектор неизвестных). Обозначим $$A = L + D + U$$, где $$L$$ — нижняя треугольная матрица с нулевыми диагональными элементами; $$D$$ — диагональная матрица; $$U$$ — верхняя треугольная матрица с нулевыми диагональными элементами.

Для решения этой системы рассмотрим итерационный процесс $${{x}^{i+1}}={{x}^{i}}-{{H}_{i}}(A{{x}^{i}}-f),$$ где $$\det ({{H}_{i}})\ne 0$$, или $${{x}^{i+1}}={{P}_{i}}{{x}^{i}}+{{d}_{i}},\qquad x(0)={{x}_{0}},$$ где $${{P}_{i}}=I-{{P}_{i}}A$$ — оператор $$i$$-го шага итерационного процесса;$${{d}_{i}}={{P}_{i}}A$$.

Итерационный процесс сходящийся, если последовательность $$\left\{ {{x}_{i}}\right\}$$ сходится к решению $${{x}^{*}}$$ при любом $$x_0$$.

Если матрица не зависит от номера итерации, итерационный процесс называется стационарным:

Необходимым и достаточным условием сходимости стационарного процесса является выполнение условия $$\rho (P)<1$$, где $$\rho (P)$$ — спектральный радиус матрицы $$P$$ (наибольшее по модулю собственное число матрицы $$P$$).

С использованием введённых обозначений метод простой итерации (метод Якоби) даётся формулой: $$P=I-{{D}^{-1}}A, \qquad D = diag({a_{ii}}), \qquad H={{D}^{-1}},$$ а метод Гаусса–Зейделя — формулой: $$P = - {(D + L)^{ - 1}}U, \qquad H={{(D+L)}^{-1}}.$$

Рассмотрим поэлементные расчетные соотношения для методов Якоби и Гаусса—Зейделя.

Все элементы главной диагонали матрицы $$I = D^{-1}A$$ равны нулю, остальные элементы равны $$-\cfrac{{a}_{ij}}{{a}_{ii}}, \,\,\, i,j =\overrightarrow {1,n}$$. Свободный член уравнения (8.4) равен $$\cfrac{{{f}_{i}}}{{{a}_{ii}}$$.

Таким образом, для метода Якоби итерационный процесс записывается в виде $$x^{k+1} = Cx^{k} +$$, где $${{C}_{ij}}=-\cfrac{{{a}_{ij}}}{{{a}_{ii}}};\ k=0,1, \dots, \,\, i,j = \overrightarrow {1,n};\ {E_i} = \cfrac{{{f_i}}}{{{a_{ii}}}}$$.

Для метода Гаусса–Зейделя $$x^{k+1}=-(D+L)^{-1}Ux^k + (D+L)^{-1}f$$, или $$(D + L){x^{k + 1}} = - U{x^k} + b$, ${x^{k + 1}} = B{x^{k + 1}} + E{x^k} +e$$, где

Рассмотрим решение конкретной системы уравнений $$Ax = b$$ методом Якоби:

Вычисляем элементы матрицы $$B$$ и вектора $$e$$:

Вычислим значения x по формуле $$x^{k+1} = Bx^k + e$$. Для решения использована следующая последовательность команд Maxima:

(%i1)	A:matrix([8,-1,2],[1,9,3],[2,-3,10])$
	b:matrix([8],[18],[5])$
	A0:matrix([A[1,1],A[1,1],A[1,1]],
		[A[2,2],A[2,2],A[2,2]],
		[A[3,3],A[3,3],A[3,3]])$
	B:-A/A0+diagmatrix(3,1)$
	e:b/matrix([A[1,1]],[A[2,2]],[A[3,3]])$ x:e$

(%i7)	xt:float(B.x+e)$
	xt:float(B.xt+e)$
	xt:float(B.xt+e)$
	xt:float(B.xt+e)$
	xt:float(B.xt+e)$
	xt:float(B.xt+e)$
	xt:float(B.xt+e)$
	x0:xt$
	xt:float(B.xt+e)$
	x1:xt$
	r:x1-x0$
	float(r); /* оценка сходимости*/
	float(A.x1-b); /* оценка невязки*/

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

Для решения системы линейных алгебраических уравнений (8.1) итерационным методом её необходимо привести к нормальному виду:

Стационарное итерационное правило получаем, если матрица $$B$$ и вектор $$\bar{g}$$ не зависят от номера итерации: $${{\bar{x}}^{k+1}}=P{{\bar{x}}^{k}}+\bar{g}$$. Нестационарное итерационное правило получаем если матрица $$B$$ или вектор $$\bar{g}$$ изменяются с ростом номера итерации: $${{\bar{x}}^{k+1}}={{P}_{k}}{{\bar{x}}^{k}}+{{\bar{g}}^{k}}$$.

Стационарное итерационное правило обычно называют методом простой итерации. Предел итерационной последовательности является точным решением системы (8.5) или (8.1).

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

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

Для определения скорости сходимости можно воспользоваться следующей теоремой: если какая-либо норма матрицы $$P$$, согласованная с данной нормой вектора, меньше единицы, то имеет место следующая оценка погрешности метода простой итерации: $$\left\| {{{\bar{x}}}^{*}}-{{{\bar{x}}}^{k}} \right\|<{{\left\| P \right\|}^{k}}\cdot \left\| {{{\bar{x}}}^{k}} \right\|+\frac{{{\left\| P \right\|}^{k}}\cdot \left\| {\bar{g}} \right\|}{1-\left\| P \right\|},$$ где $${\bar{x}}^{*}$$ — точное решение системы (8.1).

Другими словами, условие сходимости выполняется, если выполняется условие доминирования диагональных элементов матрицы исходной системы $$A$$ по строкам или столбцам: $$\sum\limits_{i,j=1}^{n}{\left| {{a}_{ij}} \right|}<\left| {{a}_{ii}} \right|$$ или $$\sum\limits_{i,j=1}^{n}{\left| {{a}_{ij}} \right|}\le \left| {{a}_{jj}} \right|$$.

В этом случае легко можно перейти от системы вида (8.1) к системе (8.5). Для этого разделим $$i$$–ое уравнение системы на $$a_{i,i}$$ и выразим $$x_i$$: $${{x}_{i}}=\frac{{{f}_{i}}}{{{a}_{ii}}}-\frac{{{a}_{11}}}{{{a}_{ii}}}{{x}_{1}} -\ldots -\frac{{{a}_{1n}}}{{{a}_{ii}}}{{x}_{n}},$$ т.е. для матрицы $$P$$ будет выполнено одно из условий сходимости, где $$P=\left[ \begin{matrix} 0 -\cfrac{{{a}_{12}}}{{{a}_{11}}} \ldots -\cfrac{{{a}_{1n}}}{{{a}_{11}}} \\ -\cfrac{{{a}_{21}}}{{{a}_{22}}} 0 \ldots -\cfrac{{{a}_{2n}}}{{{a}_{22}}} \\ \ldots \ldots \ldots \ldots \\ -\cfrac{{{a}_{n1}}}{{{a}_{nn}}} -\cfrac{{{a}_{n2}}}{{{a}_{nn}}} \ldots 0 \\ \end{matrix} \right].$$

Пример реализации метода простой итерации средствами Maxima с печатью промежуточных результатов представлен в скрипте ниже:

(%i1)	iterpr(a0,b0,x,n,eps):=block([a,b,x0,i,j,s,sum,p],
	a:copymatrix (a0), b:copymatrix(b0), x0:copymatrix(x),
	sum:1, p:0,
	while sum>eps do (
		sum:0, p:p+1, print("p= ",p," x= ",float(x)),
		for i:1 thru n do (
			s:b[i,1],
			for j:1 thru n do (s:s-a[i,j]*x0[j,1]),
			s:s/a[i,i], x[i,1]:x0[i,1]+s, sum:sum+abs(s)
		),
		x0:copymatrix(x)
	),
	float(x))$

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

В методе Зейделя система (8.1) также приводится к системе (8.5). Но при вычислении последующей компоненты вектора используются уже вычисленные компоненты этого вектора.

Итерационная формула метода в скалярной форме записывается следующим образом:

Установим связь между методом Зейделя и методом простой итерации. Для этого матрицу $$B$$ представим в виде суммы двух матриц: $$P = H + F$$, где

Итерационная формула метода Зейделя в матричной форме записывается в виде: $$\begin{align*} {{\bar{x}}^{(k+1)}}=H{{\bar{x}}^{(k+1)}}+F{{\bar{x}}^{(k)}}+\bar{g}, \qquad \text{или} \\ \left( E-H \right){{\bar{x}}^{(k+1)}}=F{{\bar{x}}^{(k)}}+\bar{g}, \qquad \text{откуда} \\ {{\bar{x}}^{(k+1)}}={{\left( E-H \right)}^{-1}}F{{\bar{x}}^{(k)}}+{{\left( E-H \right)}^{-1}}\bar{g}, \end{align*}$$ т.е. метод Зейделя эквивалентен методу простой итерации с матрицей $$(E - H)^{-1}F$$.

Исходя из полученной аналогии методов Зейделя и простой итерации, можно сформулировать следующий признак сходимости метода Зейделя: для того чтобы метод Зейделя сходился, необходимо и достаточно, чтобы все собственные значения матрицы $$(E - H)^{-1}F$$ по модулю были меньше единицы.

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

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

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

Пример реализации метода Зейделя:

(%i1)	seidel(a0,b0,x,n,eps):=block([a,b,i,j,s,sum,p],
	a:copymatrix (a0), b:copymatrix(b0),
	sum:1, p:0,
	while sum>eps do (
		sum:0, p:p+1, print("p= ",p),
		for i:1 thru n do
		(
			s:b[i,1],
			for j:1 thru n do (s:s-a[i,j]*x[j,1]),
			s:s/a[i,i], x[i,1]:x[i,1]+s, sum:sum+abs(s)
		)
	),
	float(x))$

Пример решения простой системы методом Зейделя:

(%i2)	aa:matrix([3,1,1],[1,3,1],[1,1,3]); bb:matrix([6],[6],[8]);
	x:matrix([3],[3],[3]); zz:seidel(aa,bb,x,3,0.0000001);

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

8.5.1 Методы решения задачи Коши

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

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

Пусть на отрезке $$x_0 \leq x \leq L$$ требуется найти решение $$y(x)$$ дифференциального уравнения $$y' = f(x,y),$$ удовлетворяющее при $$x = x_0$$ начальному условию $$y(x_0) = y_0$$.

Будем считать, что условия существования и единственности решения поставленной задачи Коши выполнены.

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

Отрезок $$[x_0,L]$$ накрывается сеткой (разбивается на интервалы) чаще всего с постоянным шагом $$h$$ ($$h = x_{n+1} - x_{n}$$), и по какому- то решающему правилу находится значение $$y_{n+1} = y(x_{n+1})$$. Таким образом, результатом решения задачи Коши численными методами является таблица, состоящую из двух векторов: $$x = (x_0 , x_1 , \dots , x_n)$$ — вектора аргументов и соответствующего ему вектора значений искомой функции $$y = ( y_0 , y_1, \dots , y_n )$$.

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

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

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

Приближенное решение $$y_m(x)$$ исходной задачи ищут в виде

Здесь $${y^{(0)}}(x_0) = y(x_0),\quad {y^{(1)}}(x_0) = y'(x_0) =f(x_0,y_0)$$, а значения $$y^{(i)}(x_0), i = 2, 3, \dots , m$$ находят по формулам, полученным последовательным дифференцированием заданного уравнения:

Для значений $$x$$, близких к $$x_0$$, метод рядов (8.7) при достаточно большом значении $$m$$ дает обычно хорошее приближение к точному решению $$y(x)$$ задачи (8.6). Однако с ростом расстояния $$|x - x_0|$$ погрешность приближения искомой функции рядом возрастает по абсолютной величине (при одном и том же количестве членов ряда), и правило (8.7) становится вовсе неприемлемым, когда x выходит из области сходимости соответствующего ряда (8.7) Тейлора.

Если в выражении (8.7) ограничиться $$m = 1$$, то для вычисления новых значений $$y(x)$$ нет необходимости пересчитывать значение производной, правда и точность решения будет невысока.

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

Рассмотрим интегрирование единичного дифференциального уравнения $$\cfrac{dy}{dx}=f(x,y)$$ на отрезке $$[x_0,x]$$ с начальным условием $$y(x_0) = y_0$$. При формальном интегрировании получим: $$\int\limits_{x_0}^{x}\frac{dy}{dx}dx=\int\limits_{x_0}^{x}f(x,y)dx$$ Процедура последовательных приближений метода Пикара реализуется согласно следующей схеме

В качестве примера рассмотрим решение уравнения $$y' = -y$$, при $$y(0) = 1, x_0 = 0$$:

(%i1)	rp:-y$ y0:1$ x0:0$ /*rp-правая часть уравнения */
(%i4)	y1:y0+integrate(subst(y0,y,rp),x,x0,x);
	/* y1 - первое приближение */

(%i5)	y2:y0+integrate(subst(y1,y,rp),x,x0,x);
	/* y2 - второе приближение */

(%i6)	y3:y0+integrate(subst(y2,y,rp),x,x0,x);

(%i7)	expand(%); /* Очевидное решение рассматриваемого ОДУ -
		экспонента y=exp(-x).
		В результате использование метода Пикара
		получаем решение в виде ряда Тейлора */

8.5.2 Метод рядов, не требующий вычисления производных правой части уравнения

Естественно поставить задачу о таком усовершенствовании приведенного выше одношагового метода, которое сохраняло бы основные его достоинства, но не было бы связано с нахождением значений производных правой части уравнения $$y_m(x_{n + 1}) \approx \sum\limits_{i = 0}^m \frac{h^i}{i!}y^{(i)}(x_n) ,$$ где $$x_{n+1} = x_n + h$$.

Чтобы выполнить это условие (последнее), производные y$$y^{(i)}(x)$$, $$i = 2, 3, \dots , m$$, входящие в правую часть уравнения (8.9), можно заменить по формулам численного дифференцирования их приближенными выражениями через значение функции $$y'$$ и учесть, что $$y'(x) = f [x, y(x)]$$.

8.5.2.1 Метод Эйлера

В случае $$m = 1$$ приближенное равенство (8.9) не требует вычисления производных правой части уравнения и позволяет с погрешностью порядка $$h^2$$ находить значение $$y(x_n + h)$$ решения этого уравнения по известному его значению $$y(x_n)$$. Соответствующее одношаговое правило можно записать в виде

Это правило (8.10) впервые было построено Эйлером и носит его имя. Иногда его называют также правилом ломаных или методом касательных. Метод Эйлера имеет относительно низкий порядок точности — $$h^2$$ на одном шаге. Практическая оценка погрешности приближенного решения может быть получена по правилу Рунге.

Пример реализации метода Эйлера средствами Maxima приведён в следующем примере:

(%i1)	euler1(rp,fun,y0,x0,xend,h):=block([OK,_x,_y,_y1,rez],
		_x:x0, _y:y0, rez:[_y], OK:-1, eps:0.1e-7,
		while OK<0 do (
			if ((_x+h>xend) or (abs(_x+h-xend)<eps))
				then (h:xend-_x,_x:xend, OK:1)
				else (_x:_x+h),
			_y1:makelist(float(_y[i]+h*subst([fun[i]=_y[i],x=_x],
				rp[i])),i,1,length(_y)),rez:append(rez,[_y1]),
			_y:_y1
		),
		rez
	)$

Правые части решаемых дифференциальных уравнений передаются в функцию $$euler1$$ в списке $$rp$$. По умолчанию предполагается, что список имён зависимых переменных — $$fun$$, имя независимой переменной — $$x$$. Начальные значения независимой и зависимых переменных — список $$y_0$$ и скалярная величина $$x_0$$, граница интервала интегрирования — величина $$xend$$, шаг интегрирования — $$h$$.

Следующий пример — обращение к функции $$euler1$$. Приведено решение системы из трёх дифференциальных уравнений на интервале [0, 1] с шагом $$h = 0.1$$:

С начальными условиями $$y(0) = 1.0; v(0) = 1.0; z(0) = 0$$, решением уравнений данной системы будут функции:

(%i2)	euler1([-2*y,-5*v,3*x],[y,v,z],[1,1,0],0,1,0.1);

Проверить решение можно сравнивая графики точного решения и множества вычисленных приближенных значений. Пример последовательности команд, позволяющих выделить отдельные компоненты решения системы ОДУ и построить график точного и приближенного решения третьего уравнения системы, представлен ниже (точные решения — списки $$yf, vf, zf$$; приближенные решения — списки $$yr, vr, zr;$$ список значений независимой переменной — $$xg$$).

(%i3)	rez:euler1([-2*y,-5*v,3*x],[y,v,z],[1,1,0],0,1.0,0.1)$
	n:length(rez)$
	yr:makelist(rez[k][1], k, 1, n)$
	vr:makelist(rez[k][2], k, 1, n)$
	zr:makelist(rez[k][3], k, 1, n)$
	xg:makelist(0.1*(k-1),k,1,n)$
	yf:makelist(exp(-2*xg[k]), k, 1, n)$
	vf:makelist(exp(-5*xg[k]), k, 1, n)$
	zf:makelist(3*xg[k]^2/2, k, 1, n)$
	plot2d ([[discrete, xg,zr],[discrete, xg,zf]],
			[style, points, lines])$

Уменьшение шага $$h$$ приводит к уменьшению погрешности решения (в данном примере — шаг 0.1).

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

Изложим идею метода на примере задачи Коши: $$y'=f(x,y); \\ {{x}_{0}}\le x\le b; \\ y({{x}_{0}})={{y}_{0}}. \\$$ Интегрируя это уравнение в пределах от $$x$$ до $$x + h$$ ($$0 < h < 1$$), получим равенство $$y(x + h) = y(x) + \int\limits_x^{x + h} {f\,[t,y(t)]dt},$$ которое посредством последнего интеграла связывает значения решения рассматриваемого уравнения в двух точках, удаленных друг от друга на расстояние шага $$h$$.

Для удобства записи выражения (8.11) используем обозначение $$\Delta y = y(x + h) - y(x)$$ и замену переменной интегрирования $$t = x + h$$. Окончательно получим: $$\Delta y=h\int\limits_{0}^{1}{f[x}+\alpha \,h,y(x+\alpha \,h)]d\alpha \,$$

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

Рассмотрим линейную комбинацию величин $$\phi_i,\ i = 0, 1, \dots, q$$, которая будет являться аналогом квадратурной суммы и позволит вычислить приближенное значение приращения $$\Delta y$$: $$\Delta y \approx \sum\limits_{i = 0}^q {{a_i}{\phi_i}} ,$$ где $$\begin{gathered} {\phi _0} = hf(x,y); \hfill \\ {\phi _1} = hf(x + {\alpha _1}h;y + {\beta _{10}}{\phi _0}); \hfill \\ {\phi _2} = hf(x + {\alpha _2}h;y + {\beta _{20}}{\phi _0} + {\beta _{21}}{\phi _1}); \hfill \\ \dots \hfill \\ \end{gathered}$$

Метод четвертого порядка для $$q = 3$$, являющийся аналогом широко известной в литературе четырехточечной квадратурной формулы "трех восьмых", имеет вид $$\Delta y \approx \frac{1}{8}({\phi _0} + 3{\phi _1} + 3{\phi _2} + {\phi _3}),$$ где $$\begin{gathered} {\phi _0} = hf({x_n},{y_n}); \hfill \\ {\phi _1} = hf({x_n} + \frac{h}{3},{y_n} + \frac{{{\phi _0}}}{3}); \hfill \\ {\phi _2} = hf({x_n} + \frac{2}{3}h,{y_n} - \frac{{{\phi _0}}}{3} - {\phi _1}); \hfill \\ {\phi _3} = hf({x_n} + h,{y_n} + {\phi _0} - {\phi _1} + {\phi _2}). \hfill \\ \end{gathered}$$

Особо широко известно другое вычислительное правило типа Рунге-Кутта четвертого порядка точности: $$\Delta y = \frac{1}{6}({\phi _0} + 2{\phi _1} + 2{\phi _2} + {\phi _3}),$$ где $$\begin{gathered} {\phi _0} = hf({x_n},{y_n}), \hfill \\ {\phi _1} = hf({x_n} + \frac{h}{2},{y_n} + \frac{{{\phi _0}}}{2}), \hfill \\ {\phi _2} = hf({x_n} + \frac{h}{2},{y_n} + \frac{{{\phi _1}}}{2}), \hfill \\ {\phi _3} = hf({x_n} + h,{y_n} + {\phi _2}). \hfill \\ \end{gathered}$$ Метод Рунге-Кутта имеет погрешность четвертого порядка ($$\sim h^4$$).

Функция Maxima, реализующая метод Рунге-Кутта 4-го порядка, приведена в следующем примере (с печатью промежуточных результатов):

(%i1)	rk4(rp,fun,y0,x0,xend,h):=block(
		[OK,n,h1,_x,_y,_k1,_k2,_k3,_k4,rez],
		_x:x0, _y:y0, rez:[_y], OK:-1, h1:h, n:length(_y),
		while OK<0 do (
			if (_x+h1>=xend) then (h1:xend-_x, OK:1),
			_k1:makelist(float(h1*subst([fun[i]=
				float(_y[i]),x=float(_x)],rp[i])),i,1,n),
			_k2:makelist(float(h1*subst([fun[i]=
				float(_y[i]+_k1[i]/2),x=float(_x+h1/2)],
					rp[i])),i,1,n),
			_k3:makelist(float(h1*subst([fun[i]=
				float(_y[i]+_k2[i]/2),x=float(_x+h1/2)],
					rp[i])),i,1,n),
			_k4:makelist(float(h1*subst([fun[i]=
				float(_y[i]+_k3[i]),x=float(_x+h1)],rp[i])),
					i,1,n),
			_y1:makelist(float(_y[i]+
				(_k1[i]+2*_k2[i]+2*_k3[i]+_k4[i])/6),i,1,n),
			rez:append(rez,[_y1]),
			print("x= ",_x, " y= ",_y),
			_x:_x+h1,
			_y:_y1
		), rez
	)$

Пример обращения к функции $$rk4$$ представлен следующей последовательностью команд (решалась та же система, что и при тестировании метода Эйлера):

(%i2)	rk4([-2*y,-5*v,3*x],[y,v,z],[1,1,1],0,1,0.1);
x= 0	y= [1,1,1]
x= 0.1	y= [0.81873333333333,0.60677083333333,1.015]
x= 0.2	y= [0.67032427111111,0.36817084418403,1.06]
x= 0.3	y= [0.54881682490104,0.22339532993458,1.135]
x= 0.4	y= [0.44933462844064,0.13554977050718,1.24]
x= 0.5	y= [0.3678852381253,0.082247647208783,1.375]
x= 0.6	y= [0.30119990729446,0.04990547343658,1.54]
x= 0.7	y= [0.24660240409888,0.030281185705008,1.735]
x= 0.8	y= [0.20190160831589,0.018373740284549,1.96]
x= 0.9	y= [0.16530357678183,0.011148649703906,2.215]
x= 1.0	y= [0.13533954843051,0.0067646754713805,2.5]