«Решение систем нелинейных уравнений»

ГОУ ОГУ
Лабораторная работа M 2.2

Тема: «Решение систем нелинейных уравнений »
Цель занятия:

После этого занятия вы должны УМЕТЬ:
1 Реализовать и применять методы решения систем нелинейных уравнений : простых итераций, Ньютона

2 Составлять схемы алгоритмов для методов : простых итераций, Ньютона
Задачи занятия:

1) Изучить основные теоретические моменты по теме.

2) Разобрать готовую программу.

3) Написать собственную программу.

4) Составить схемы алгоритмов для методов
Время: 4 часа
Порядок выполнения лабораторной работы:
1 Внимательно изучить теоретический материал.

2 В среде Турбо Паскаль проработать готовую программу по методу простых итераций (приложение А).

3 Составить блок схемы по готовой программе.

4 Ответить на контрольные вопросы.

5 Написать собственную программу по методам решения систем нелинейных уравнений и систем нелинейных уравнений и составить блок-схему для нее (вариант задания см. приложение Б).

6 Составить отчет по лабораторной работе.

7 Защитить лабораторную работу по пунктам 2-6.
Содержание отчета по лабораторной работе:

1 Тема, цель

2 Постановка задачи

3 Теоретические предпосылки (ответить на контрольные вопросы)

4 Приложение А: Текст программы

5 Приложение В: Блок-схемы

Теоретические предпосылки

1 Решение систем нелинейных уравнений
Рассмотрим систему нелинейных уравнений с m неизвестными вида

.

В отличие от систем линейных уравнений здесь использование прямых методов исключено и решение находится с использованием итерационных методов, т.е. находится приближенное решение
x^* = (x₁^*, ... , x_m^*),
удовлетворяющее при заданном  > 0 условию .

. Введем матрицу Якоби
.
Как и в случае решения одного уравнения начинаем с этапа локализации решения (отделения корней).

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

Дана система нелинейных уравнений
(1)

или в векторной форме

x = (x_{1 ,} … x_m ) ;

f = (f₁, …f_m ) ^T

f (x) = 0
Алгоритм решения систем
1. Исходную систему приводят к виду

в векторной форме

x = ( x ) (2 )

2. Определяют начальное приближение

x^{( 0 )} = (x₁ ⁽⁰⁾ , … , x_m^{( 0 )} )

3. 
   Проверяют условие сходимости метода ║φ΄(x ⁽⁰⁾ )║< 1
   

4. 
   Если условие сходимости выполняется, то каждое последующее приближение получают путем подстановки предыдущего в систему (2)
   

5. 
   Условие окончания счета
   

║x ^{( к + 1)} – x ^{( к)}║<
3 Метод Ньютона
Дана система нелинейных уравнений (1)

f (x) = 0

x = (x₁, x₂, …x_m )

f = (f₁, f₂, … f_m ) ^T
Алгоритм метода Ньютона

1. 
   Задают начальное приближение
   

x ⁽⁰⁾ = (x₁ ⁽⁰⁾, x ₂⁽⁰⁾,…x_m⁽⁰⁾)

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

∆ x ^{( к+1)} :

f ΄ (x ^{( к )} ) ∆ x ^{(к + 1 )} = -f ( x ^{( к)} ) и

полагая затем

x ^{( к+1)} = x ^{( к)} + ∆ x ^{( к+1)} ,

где

f ΄(x) =

3. 
   Критерий остановки
   

║x ^{( к+1)} – x ^{( к )} ║ <

Контрольные вопросы

1. 
   Что называется нормой матрицы? Напишите формулы для ее расчета
   
2. 
   Сформулируйте достаточные условия сходимости метода простых итераций
   
3. 
   Сформулируйте достаточные условия сходимости метода Ньютона
   
4. 
   Какова оценка погрешности и скорость сходимости метода простых итераций
   
5. 
   Какова оценка погрешности и скорость сходимости метода Ньютона
   
6. 
   В чем состоит метод простой итерации для решения системы нелинейных уравнений?
   
7. 
   В чем метод Ньютона превосходит метод простой итерации?
   
8. 
   В чем состоят основные трудности, связанные с практическим применением метода Ньютона?
   
9. 
   Какое условие окончания итераций применяется на практике?
   

ПРИЛОЖЕНИЕ А

(cos(x-1)+y=0.5)&(x-cos(y)=3)"

3.3

1.2

1 Цель занятия

1 Реализовать и применять методы : простых итераций, Ньютона.

2 Составлять схемы алгоритмов для методов: простых итераций, Ньютона.
2 Постановка задачи

Система уравнений:



Якобиан:


#include

#include

struct _Result

{

double x;

double y;

};

//---------------------------------------------------------------------------
//-----------------------------------------------------------------------------
//-----------------------------------------------------------------------------

double Function1(double x, double y)

{

return (cos(x-1)+y-0.5);

}

//---------------------------------------------------------------------------

double Function2(double x, double y)

{

return (x-cos(y)-3);

}

//**********

//---------------------------------------------------------------------------

#pragma package(smart_init)

void MakeMatrixOfSubstitutionFuntions(double xy[2],double MatrixOfSubstitutionFuntions[2])

{

MatrixOfSubstitutionFuntions[0]=3+cos(xy[1]);

MatrixOfSubstitutionFuntions[1]=0.5-cos(xy[0]-1);

}

//---------------------------------------------------------------------------

void MakeMatrixOfFuntions(double xy[2],double MatrixOfFuntions[2])

{

MatrixOfFuntions[0]=Function1(xy[0],xy[1]);

MatrixOfFuntions[1]=Function2(xy[0],xy[1]);

}

//----------------------------------------------------------------------------

void MakeJacobian(double xy[2],double Jacobian[4])

{

Jacobian[0]=-sin(xy[0]-1);

Jacobian[1]=1;

Jacobian[2]=1;

Jacobian[3]=sin(xy[1]);

}

//----------------------------------------------------------------------------

double DetMatrix(double Matrix[4])

{

return Matrix[0]*Matrix[3]-Matrix[1]*Matrix[2];

}

//----------------------------------------------------------------------------

void InverseMatrix(double SourseMatrix[4],double InverseMatrix[4])

{

if (DetMatrix(SourseMatrix)==0)

{

cout<<"Singular matrix!"<

}

else

{

double InverseDeterminat=1/DetMatrix(SourseMatrix);
InverseMatrix[0]=InverseDeterminat*SourseMatrix[3];

InverseMatrix[1]=InverseDeterminat*SourseMatrix[2]*(-1);

InverseMatrix[2]=InverseDeterminat*SourseMatrix[1]*(-1);

InverseMatrix[3]=InverseDeterminat*SourseMatrix[0];

}
}

//----------------------------------------------------------------------------

void ProductMatrixes(double MatrixA[4], double MatrixB[2],double MatrixResult[2])

{

MatrixResult[0]=MatrixA[0]*MatrixB[0]+MatrixA[1]*MatrixB[1];

MatrixResult[1]=MatrixA[2]*MatrixB[0]+MatrixA[3]*MatrixB[1];

}

//----------------------------------------------------------------------------

void DifferenceMatrixes(double MatrixA[2], double MatrixB[2],double MatrixResult[2])

{

MatrixResult[0]=MatrixA[0]-MatrixB[0];

MatrixResult[1]=MatrixA[1]-MatrixB[1];

}

//----------------------------------------------------------------------------

double NormaOfDifference(double Vector1[2],double Vector2[2])

{

double Norma=0;
for (int i=0;i<=1;i++)

{

Norma+=fabs(Vector1[i]-Vector2[i]);

}

return Norma;

}

//----------------------------------------------------------------------------

void CopyMatrix(double SourseMatrix[2], double DestMatrix[2])

{

DestMatrix[0]=SourseMatrix[0];

DestMatrix[1]=SourseMatrix[1];

}

//----------------------------------------------------------------------------

int NewtonsMetod(double Accuracy,_Result InitialApproximation, _Result &Result;)

{

double xy[2];

double xy_previous[2];

double Jacobian[4];

double MatrixOfFuntions[2];

double InverseJacobian[4];

double MatrixResult[2];

int CountOfIterations=1;
xy_previous[0]=InitialApproximation.x;

xy_previous[1]=InitialApproximation.y;
MakeJacobian(xy_previous,Jacobian);

if (DetMatrix(Jacobian)==0)

{

cout<<"Singular Jacobian!"<

return 0;

}
MakeMatrixOfFuntions(xy_previous,MatrixOfFuntions);

InverseMatrix(Jacobian,InverseJacobian);

ProductMatrixes(InverseJacobian,MatrixOfFuntions,MatrixResult);

DifferenceMatrixes(xy_previous,MatrixResult,xy);
while(NormaOfDifference(xy,xy_previous)>=Accuracy)

{

CopyMatrix(xy,xy_previous);

MakeJacobian(xy_previous,Jacobian);

if (DetMatrix(Jacobian)==0)

{

cout<<"Singular Jacobian! Accuracy isn\'t reached!"<

break;

}
MakeMatrixOfFuntions(xy_previous,MatrixOfFuntions);

InverseMatrix(Jacobian,InverseJacobian);

ProductMatrixes(InverseJacobian,MatrixOfFuntions,MatrixResult);

DifferenceMatrixes(xy_previous,MatrixResult,xy);

CountOfIterations++;

}
Result.x=xy[0];

Result.y=xy[1];

return CountOfIterations;

}

//---------------------------------------------------------------------------

int MetodOfIterations(double Accuracy,_Result InitialApproximation, _Result &Result;)

{

double xy[2];

double xy_previous[2];

double MatrixOfSubstitutionFuntions[2];

int CountOfIterations=0;
xy_previous[0]=InitialApproximation.x;

xy_previous[1]=InitialApproximation.y;
MakeMatrixOfSubstitutionFuntions(xy_previous,xy);
while(NormaOfDifference(xy,xy_previous)>=Accuracy)

{

CopyMatrix(xy,xy_previous);

MakeMatrixOfSubstitutionFuntions(xy_previous,xy);

CountOfIterations++;

}
Result.x=xy[0];

Result.y=xy[1];
return CountOfIterations;

}

//-----------------------------------------------------------------------------
#pragma hdrstop

#pragma argsused
#include "Unit2.cpp"

#include

#include

//**********

int main(int argc, char* argv[])

{
double Accuracy;

int LevelOfRounding;

_Result Result;

_Result InitialApproximation;
cout << "Equations set: (cos(x-1)+y=0.5)&(x-cos(y)=3)"<

cout << "Enter accuracy: ";

cin >> Accuracy;
InitialApproximation.x=3.3;

InitialApproximation.y=1.2;
cout<<"Newton\'s metod:"<

cout<<"Count of iterations: "<

cout<<"result x: "<

cout<<"Metod of iterations:"<

cout<<"Count of iterations: "<

cout<<"result x: "<
cin>>"\n";

return 0;

}

//---------------------------------------------------------------------------

Метод простых итераций
Порядок выполнения:

1. 
   Найти нулевое приближение решения;
   
2. 
   Преобразовать систему f(x) = 0 к виду х = Ф(х);
   
3. 
   Проверить условие сходимости метода.
   

Схема алгоритма приведена на рисунке 1.

Пример. Решить методом простых итераций систему

В качестве нулевого приближения выберем точку х=1, у = 2.2, z = 2. Преобразуем систему к виду


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

uses Crt;
type mas=array [0..1] of comp;
var

x10,x20,e,max: real;

fii,x:mas;

{____________________________________}

function f(x1,x2: real): real;

begin

f:=x1*ln(x2)-x2*ln(x1);

end;

{____________________________________}

function fi(x1,x2: real;var fii:array of comp): real;

begin

fii[0]:=exp(1/3*ln(8*x1*x2-(x2*x2*x2)));

fii[1]:=x2+x2/ln(x2)-x1/ln(x1);

end;

{____________________________________}

function iter(x10,x20: real;var x:array of comp): real;

var x1,x2:real;

n:integer;

begin

n:=0;
repeat

x1:=x10;

x2:=x20;
n:=n+1;
fi(x10,x20,fii);

x10:=fii[0];

x20:=fii[1];

max:=abs(x10-x1);

if abs(x20-x2)>max then max:=abs(x20-x2);

until max

writeln('e = ',e:20:6);

writeln(' KOLICHESTWO ITERAZIJ = ',n);

x[0]:=x10;

x[1]:=x20;

writeln(' x1 = ',x[0]:20:4);

writeln(' x2 = ',x[1]:20:4);

end;

{___________________________________}

begin

ClrScr;
e:=1e-3;

x10:=3.8;

x20:=2;

iter(x10,x20,x);
writeln(' * * * PROWERKA * * * ');

e:=f(x[0],x[1]);

writeln(' f2 = ',e:20:4);

readln

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


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

Программу можно использовать для решения систем не выше десятого порядка.

Входные параметры: n — число уравнений системы (совпадает с числом неизвестных), n£10; х—массив из n действительных чисел, содержащий начальное приближение решения; f— имя внешней процедуры f(n, х, у), вычисляющей по заданным значениям х, расположенным в элементах массива х, текущие значения функции f и размещающей их в элементах массива у; g — имя внешней процедуры g(n, x, d), вычисляющей по заданным значениям х из массива х элементы матрицы ,которая размещается в массиве d размерности n x n; eps — значение условия окончания итерационного процесса.

Выходные параметры: х — массив из n действительных чисел (он же входной) содержит при выходе из подпрограммы приближенное значение решения; k—количество итераций.

Порядок выполнения:

1. 
   Определить из геометрических соображений начальное приближение решения.
   
2. 
   Составить подпрограммы f и g для вычисления значений f и.
   
3. 
   Составить головную программу, содержащую описание массива х, имен f и g, присвоение фактических значений параметрам n, х, eps, обращение к ТNUTON и печать результатов.
   

4. Провести вычисления на ЭВМ.
Пример. Решить систему уравнений


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

X(1)= 1.1132 Х(2)= 2.3718 Х(3)= 2.1365

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

ПРИЛОЖЕНИЕ Б

ПРИЛОЖЕНИЕ Б₁

Варианты заданий

Таблица Б₁

1			16
2			17
3			18	sin(y+l)-x=l,2 2y+cos x = 2.
4	sin у+2х=2 cos(x-1)+y = 0,7		19	sin(y+0,5)-x=l cos(x-2)+y =0
5	sin(y-l)+x=l,3 y-sin(x+l)=0,8		20	cos(y+0,5)-x = 2 sin x —2y= 1
6	sin(x+l)-y=l 2x+cos y = 2		21	sin x+2y= 1,6 cos(y-1)+x=l
7	sin(x+0,5)-y=1,2 cos(y-2)+x = 0		22	cos(x-l)+y = 0,5 x –cos y = 3
8	cos x+y= 1,5 2x-sin(y-0,5)=l		23	cos (x+0,5)+у = 0,8 sin y — 2x=l,6
9	2y-cos(x+l)=0 x+sin y= -0,4		24
10	cos(y-l)+x=0,5 y-cos x = 3		25	cos y+x=l,5 2y-sin(x-0,5)=l
11	cos(y+0,5) + x=0,8 sin x-2y=l,6		26	2x-cos(y+l) = 0 y+sin x= -0,4
12	sin(y + 2)-x=l,5 y+cos(x-2)=0,5.		27	cos(x-l)+y=0,8 x—cos y=2
13	cos x+y= 1,2 2x-sin(y-0,5) = 2.		28	cos(x+0,5)+y=l siny—2x=2
14	sin(x-1)+y= 1,5 х - sin(y+ 1)= 1		29	sin(у+1)-х=1 2y+cos x = 2
15	cos(y-l)+x=0,8 y-cos x=2.		30	cos(x-l)+y=l sin y+2x= 1,6