1.4 Решение систем нелинейных уравнений. Метод Ньютона
Для системы второго порядка
f(x,y) = 0,
g(x,y) = 0
последовательные приближения по методу Ньютона вычисляются по формулам
,
,
где
,
,
.
Метод Ньютона сходится, если начальное приближение (x0, y0) выбрано удачно и матрица Якоби невырожденная, причём сходимость квадратичная. На практике итерации обычно заканчивают, когда одновременно достаточно малы значения |f(xn,yn)| и |g(xn,yn)| или разности |xn+1–xn| и |yn+1–yn|. Для выбора начального приближения применяют графический метод, метод проб и т.д.
Пример. Решить систему уравнений
(1.12)
Решение в системе Maple:
1) Очистим оперативную память. Нам в дальнейшем потребуются функции det и implicitplot, которых нет в ядре системы, поэтому здесь подключим пакеты линейной алгебры linalg и графических построений plots.
> restart;
> with(linalg): with(plots):
2) Зададим функции f(x,y) и g(x,y).
> f:=(x,y)->evalf(x^7-5*x^2*y^4+1510);
> g:=(x,y)->evalf(y^5-3*x^4*y-105);
3) Для отделения корней системы построим графики неявных функций f(x,y)=0 и g(x,y)=0.
> implicitplot({f(x,y)=0,g(x,y)=0},x=-5..5,y=-5..5);
Рис. 1.4.
Видим пять пересечений графиков функций, значит, система имеет как минимум пять решений.
4) Определим якобиан
> J:=[[diff(f(x,y),x),diff(f(x,y),y)], [diff(g(x,y),x),diff(g(x,y),y)]];
5) Введём функцию Jn(x,y), являющуюся определителем матрицы J.
> Jn:=unapply(det(J),x,y);
6) Определим величины An и Bn в виде соответствующих определителей, при этом для исключения символьных преобразований результат оформляем через evalf
> An:=evalf(det([[f(x,y),diff(f(x,y),y)], [g(x,y),diff(g(x,y),y)]]) );
> Bn:=evalf(det([[diff(f(x,y),x),f(x,y)], [diff(g(x,y),x),g(x,y)]]) );
7) Присвоим начальные значения переменным. x0 и y0 определяет точку начального приближения (определяется графически).
> x0:=-3.; y0:=-1.; eps:=0.001: n:=0:
8) Организуем цикл вычислений по методу Ньютона
> while abs(evalf(f(x0,y0)))+abs(evalf(g(x0,y0)))>eps do
if abs(evalf(Jn(x0,y0)))<10^(-4) then
print(`Якобиан = 0`); break fi; # Если якобиан = 0, выйти из программы
x0:=x0-subs(x=x0,y=y0,An)/Jn(x0,y0);
y0:=y0-subs(x=x0,y=y0,Bn)/Jn(x0,y0);
if n=50 then break fi; # Защита от зацикливания
n:=n+1;
od:
9) Вывод результата: (x0,y0) – решение системы, n – число итераций
> evalf({x0,y0});n;
{–2.844483289, –.5348543094}
4
10) Для вычисления остальных четырёх корней системы нужно выполнить шаги 7–9, предварительно изменив в п. 7 начальное приближение. В частности, установив
> x0:=-2.; y0:=-3.;
получим второе решение системы:
{–2.573256593, –2.304767678}
4
Третье, четвертое и пятое решение получим с начальными приближениями а) x0:=2., y0:=–3., б) x0:=-2., y0:=3., в) x0:=2., y0:=3. соответственно. Эти значения следует последовательно подставить в п. 7 и выполнить все последующие команды программы.
В завершение сведем все найденные решения в ответе.
Ответ: (–2,844; –0,535), (–2,573, –2,305), (2,430; –2,873), (–1,922, 2,957), (1,994; 3,012).
Из других методов решения систем нелинейных уравнений имеют применение метод простой итерации, метод наискорейшего спуска, градиентные методы.