К содержанию

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

> restart;

> epsilon:=0.00001; Задаем погрешность вычисления.

> with(linalg): Активизация пакета линейной алгебры.

Warning, the protected names norm and trace have been redefined and unprotected
 

> eqn1:=tan(x-y)=x*y; Задание системы уравнений и задание функций.

> eqn2:=x^2+2*y^2=1;

> f1:=(x,y)->sin(x-y)/cos(x-y)-x*y;

> f2:=(x,y)->x^2+2*y^2-1;

Задается нулевое приблежение

> x0:= Pi/3;

> y0:= 1;

> X0:=matrix(2,1,[x0,y0]);

>

> X:=matrix(2,1,[x,y]);

> F:= matrix(2,1,[f1(x,y),f2(x,y)]);

Задаем функции, являющиеся элемнтами матрицы Якоби:

> J11:=(x,y)->1/(cos(x-y))^2-y;

> J12:=(x,y)->-1/(cos(x-y)^2)-x;

> J21:=(x,y)->2*x;

> J22:=(x,y)->4*y;

> J0:=matrix(2,2,[J11(x0,y0),J12(x0,y0),J21(x0,y0),J22(x0,y0)]):evalf(%);

Достаточное условие сходимости:

> is(norm(J0)<>0);

> F0:= matrix(2,1,[f1(x0,y0),f2(x0,y0)]):evalf(%);

> X1=evalm(X0-1/J0*F0): evalf(%);

> k:=0; r:= 1000*epsilon;

Цикл для выполнения метода:

> while r>epsilon do X0:= matrix(2,1,[x0,y0]);  J0:=evalf(matrix(2,2,[J11(x0,y0),J12(x0,y0),J21(x0,y0),J22(x0,y0)])); F0:= evalf(matrix(2,1,[f1(x0,y0),f2(x0,y0)])); X1:=evalf(evalm(X0-1/J0*F0)); r:= evalf(abs(norm(evalm(X1-X0)))); k:=k+1:print('_____________СЛЕДУЮЩАЯ_ИТЕРАЦИЯ______________'); x0:=X1[1,1]: y0:=X1[2,1]: end;

...

где r - погрешность вычисления на каждом шаге итерации (определяет остановку цикла)!

истинные корни:

> solve({eqn1,eqn2},{x,y});

Проверка:

> is(abs(.7775316909-x0)<epsilon);

> is(abs(.4446596842-y0)<epsilon);

Наверх

К содержанию