Итерационные методы решения систем нелинейных уравнений (стр. 7 из 8)

z0:=StrToFloat(Edit3.Text);

eps:=StrToFloat(Edit20.Text);

i:=1;

x:=0.1-x0*x0+2*y0*z0;

y:=-0.2+y0*y0-3*x0*z0;

z:=0.3-z0*z0-2*x0*y0;

repeat

i:=i+1;

x0:=x;

y0:=y;

z0:=z;

x:=0.1-x0*x0+2*y*z;

y:=-0.2+y0*y0-3*x0*z0;

z:=0.3-z0*z0-2*x0*y0;

until ((abs(x-x0)<eps)and(abs(y-y0)<eps)and(abs(z-z0)<eps));

Edit8.Text:=FloatToStr(x);

Edit9.Text:=FloatToStr(y);

Edit10.Text:=FloatToStr(z);

Edit11.Text:=IntToStr(i);

end;

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

Строим матрицу Якоби:

> restart;

> with(LinearAlgebra):

> f1:=0.1-x0^2+2*y0*z0-x0;

> f2:=-0.2+y0^2-3*x0*z0-y0;

> f3:=0.3-z0^2-2*x0*y0-z0;

> f1x:=diff(f1,x0);

> f1y:=diff(f1,y0);

> f1z:=diff(f1,z0);

> f2x:=diff(f2,x0);

> f2y:=diff(f2,y0);

> f2z:=diff(f2,z0);

> f3x:=diff(f3,x0);

> f3y:=diff(f3,y0);

> f3z:=diff(f3,z0);

> A:=<<f1x|f1y|f1z>,<f2x|f2y|f2z>,<f3x|f3y|f3z>>;

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

> x0:=0;y0:=0;z0:=0;

> A:=A;

> A1:=A^(-1);

> f:=<f1,f2,f3>;

> X0:=<x0,y0,z0>:

> X:=Add(X0,(Multiply(A1,f)),1,-1);

> X0:=X;

> x0:=X[1];y0:=X[2];z0:=X[3];

> A:=<<f1x|f1y|f1z>,<f2x|f2y|f2z>,<f3x|f3y|f3z>>;

> A1:=A^(-1);

> f:=<f1,f2,f3>;

> X:=Add(X0,(Multiply(A1,f)),1,-1);

> i:=2:

> while (Norm(f))>0.0001 do

X0:=X;

x0:=X[1];y0:=X[2];z0:=X[3];

A:=<<f1x|f1y|f1z>,<f2x|f2y|f2z>,<f3x|f3y|f3z>>;

A1:=A^(-1);

f:=<f1,f2,f3>;

X:=Add(X0,(Multiply(A1,f)),1,-1);

i:=i+1;

end do:

> X:=X;

Получилиответ:

Количество итераций:

Текстпрограммы

procedure TForm1.Button4Click(Sender: TObject);

type mas=array[1..3]of real;

matr=array[1..3,1..3]of real;

var a,a1:matr;

i,j:integer;

x,y,z,eps:real;

f,xk,xkk,temp:mas;

procedure jacobi(x,y,z:real; var a:matr);

begin

a[1,1]:=-2*x-1;

a[1,2]:=2*z;

a[1,3]:=2*y;

a[2,1]:=-3*z;

a[2,2]:=2*y-1;

a[2,3]:=-3*x;

a[3,1]:=-2*y;

a[3,2]:=-2*x;

a[3,3]:=-2*z-1;

end;

procedure inverse(a:matr;var a1:matr);

var i,j:integer;

det:real;

s:matr;

begin

det:=a[1,1]*a[2,2]*a[3,3]+a[2,1]*a[3,2]*a[1,3]+a[1,2]*a[2,3]*a[3,1]-a[3,1]*a[2,2]*a[1,3]-a[3,2]*a[2,3]*a[1,1]-a[2,1]*a[1,2]*a[3,3];

s[1,1]:=a[2,2]*a[3,3]-a[2,3]*a[3,2];

s[1,2]:=-a[1,2]*a[3,3]+a[1,3]*a[3,2];

s[1,3]:=a[1,2]*a[2,3]-a[1,3]*a[2,2];

s[2,1]:=-a[2,1]*a[3,3]+a[2,3]*a[3,1];

s[2,2]:=a[1,1]*a[3,3]-a[1,3]*a[3,1];

s[2,3]:=-a[1,1]*a[2,3]+a[1,3]*a[2,1];

s[3,1]:=a[2,1]*a[3,2]-a[2,2]*a[3,1];

s[3,2]:=-a[1,1]*a[3,2]+a[1,2]*a[3,1];

s[3,3]:=a[1,1]*a[2,2]-a[1,2]*a[2,1];

for i:=1 to 3 do

for j:=1 to 3 do

a1[i,j]:=(1/det)*s[i,j];

end;

procedure fx(x,y,z:real; var f:mas);

begin

f[1]:=-x-x*x+2*y*z+0.1;

f[2]:=-y+y*y-3*x*z-0.2;

f[3]:=-z-z*z-2*x*y+0.3;

end;

procedure minus(x,y:mas; var z:mas);

var i:integer;

begin

for i:=1 to 3 do

z[i]:=x[i]-y[i];

end;

function max(f:mas):real;

var p:real;

i:integer;

begin

p:=0;

for i:=1 to 3 do

if abs(f[i])>p then p:=abs(f[i]);

max:=p;

end;

procedure mult(a:matr;b:mas;var c:mas);

begin

c[1]:=a[1,1]*b[1]+a[1,2]*b[2]+a[1,3]*b[3];

c[2]:=a[2,1]*b[1]+a[2,2]*b[2]+a[2,3]*b[3];

c[3]:=a[3,1]*b[1]+a[3,2]*b[2]+a[3,3]*b[3];

end;

begin

xk[1]:=StrToFloat(Edit1.Text);

xk[2]:=StrToFloat(Edit2.Text);

xk[3]:=StrToFloat(Edit3.Text);

eps:=StrToFloat(Edit20.Text);

i:=0;

repeat

fx(xk[1],xk[2],xk[3],f);

jacobi(xk[1],xk[2],xk[3],a);

inverse(a,a1);

mult(a1,f,temp);

minus(xk,temp,xk);

i:=i+1;

until max(f)<eps;

Edit12.Text:=FloatToStr(xk[1]);

Edit13.Text:=FloatToStr(xk[2]);

Edit14.Text:=FloatToStr(xk[3]);

Edit15.Text:=IntToStr(i);

end;

3.4 Модифицированный метод Ньютона

Аналогично методу Ньютона построим матрицу Якоби для данной системы уравнений, выберем начальное приближение заведомо близко к решению, построим последовательность:

> with(LinearAlgebra):

> f1x:=diff(f1,x0);

> f1y:=diff(f1,y0);

> f1z:=diff(f1,z0);

> f2x:=diff(f2,x0);

> f2y:=diff(f2,y0);

> f2z:=diff(f2,z0);

> f3x:=diff(f3,x0);

> f3y:=diff(f3,y0);

> f3z:=diff(f3,z0);

> A:=<<f1x|f1y|f1z>,<f2x|f2y|f2z>,<f3x|f3y|f3z>>;

> x0:=0;y0:=0;z0:=0;

> A:=A;

> A1:=A^(-1);

> f:=<f1,f2,f3>;

> X0:=<x0,y0,z0>;

> X:=Add(X0,(Multiply(A1,f)),1,-1);

> X0:=X;

> x0:=X[1];y0:=X[2];z0:=X[3];

> f:=<f1,f2,f3>;

> i:=1;

> while (Norm(f))>0.0001 do

X0:=X;

x0:=X[1];y0:=X[2];z0:=X[3];

f:=<f1,f2,f3>;

X:=Add(X0,(Multiply(A1,f)),1,-1);

i:=i+1;

end do;

Получилиответ:

Количество итераций:

Текстпрограммы

procedure TForm1.Button2Click(Sender: TObject);

type mas=array[1..3]of real;

matr=array[1..3,1..3]of real;

var x,y,z,ex,ey,ez,eps,b,c,d:real;

r,xk,f,temp:mas;

a,h,w,a1:matr;

i,kk: integer;

procedure jacobi(x,y,z:real; var a:matr);

procedure inverse(a:matr;var a1:matr);

procedure fx(x,y,z:real; var f:mas);

procedure minus(x,y:mas; var z:mas);

function max(f:mas):real;

procedure mult(a:matr;b:mas;var c:mas);

// все процедуры полностью совпадают с описанными выше реализации метода Ньютона

begin

xk[1]:=StrToFloat(Edit1.Text);

xk[2]:=StrToFloat(Edit2.Text);

xk[3]:=StrToFloat(Edit3.Text);

eps:=StrToFloat(Edit20.Text);

i:=0;

jacobi(xk[1],xk[2],xk[3],a);

inverse(a,a1);

repeat

fx(xk[1],xk[2],xk[3],f);

mult(a1,f,temp);

minus(xk,temp,xk);

i:=i+1;

until max(f)<eps;

Edit16.Text:=FloatToStr(xk[1]);

Edit17.Text:=FloatToStr(xk[2]);

Edit18.Text:=FloatToStr(xk[3]);

Edit19.Text:=IntToStr(i);

Выводы

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

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

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

На нашем примере можно в этом убедиться.

Метод градиентного спуска при начальном приближении даже равном

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

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

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

. При начальном приближении итерационная последовательность ещё сходится к решению, если взять начальные условия дальше от корня – метод не сходится ни к указанному выше, какому другому решению системы. Исходя из имеющихся данных про точное решение системы нелинейных уравнений, мы строим последовательность. Эта последовательность сходится только с начальным приближением лежащим в окрестности выбранного корня, и только к нему. Несмотря на возможную близость начального приближения к какому-то другому решению, она к нему не сойдётся, в отличие от других итерационных методов. Это говорит о том, что для каждого из множества решений системы нужно строить своё отображение, удовлетворяющее условиям сходимости. В этом проявляется недостаток метода простых итераций. Но если сжимающее отображение построено правильно, то преимущество метода состоит в простоте вычислений.

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