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

Лабораторная работа 4

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

 

 

Вариант 2.

Решить систему уравнений  методом Ньютона и его упрощением с точностью ε=10^-6:

- исследовать устойчивость  и сходимость методов

- сравнить методы между  собой

- начальное приближение  найти графически.

Система:

 

 

 

 

 

 

 

 

 

 

 

Теория.

Системы нелинейных уравнений

                            

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

                                     ,

где – решение системы линейных уравнений

                   

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

   Итерации методом  Ньютона продолжают до тех  пор, пока не будет получено  , где – заданная точность решения.

   Нулевое приближение  в случае двух переменных можно найти графически: построить на плоскости кривые и и найти точки их пересечения. Для трёх и более переменных удовлетворительных способов подбора нулевых приближений нет.

 

Задание 1 .Решить систему уравнений методом Ньютона с точностью ε=0,000001:

Система:

Перепишем систему в виде:

X=sin(y+1)-1,2

Y=(3+cos(x))/2

Графически находим  начальное приближение    x⁽⁰⁾=-1; y⁽⁰⁾=2

 

 

Метод Ньютона  запишем в координатной форме:

                  ,

                  .

 

 

Перепишем систему в следующем  виде:

f₁(x,y)= sin(y+1)-x -1,2=0

f₂(x,y)=2y-cos x-3=0

Найдём частные производные:

(f₁)^’ _x (x,y)= -1;

(f₂ )^’_x(x,y)=sin(x);

-------------------------------------

(f₁)^’ _y(x,y)= cos(y+1)-;

(f₂ )^’_y(x,y)=2;

Матрица Якоби:

-1         cos(y+1)

sin(x)       2

 

Окончательный (рабочий) вид системы:

-1∆x^(k)+cos(y+1)∆y^(k)= -f₁(x^(k),y^(k)),

sin(x^(k))∆x^(k)+2∆y^(k)= -f₂(x^(k),y^(k)

Теперь можно  начинать итерации. Для первой итерации требуются значения и   (x⁽⁰⁾=-1; y⁽⁰⁾=2) .Их находим, решая систему:

-1∆x⁽⁰⁾-0,9900∆y⁽⁰⁾= 0.0589

-0,8415*∆x⁽⁰⁾+2∆y⁽⁰⁾= -0,4597

∆x⁽⁰⁾=0.1191

∆y⁽⁰⁾=-0,1798

--------------------------------------------------------

x⁽¹⁾=x⁽⁰⁾+ ∆x⁽⁰⁾=-1+0.1191=-0,8809;

y⁽¹⁾=1,8202

∆x⁽¹⁾=-0.0010

∆y⁽¹⁾=-0.0024

---------------------------------------------------------------

x⁽²⁾ =-0,8819

y⁽²⁾ =1,8178

Δy⁽²⁾= -0,0000004

Δx⁽²⁾= -0,0000006

-------------------------------------

x⁽³⁾ =-0,8819

y⁽³⁾ =1,8178

 

Неравенство выполняется, поэтому итерации прекращаем.

ОТВЕТ: x⁽³⁾ =-0,8819

                  y⁽³⁾ =1,8178

EXCEL:

MATLAB:

 

x1=-2:0.1:2;

y1=sin(x1-0.6)-1.6;

y2=-3:0.1:3;

x2=(0.9+cos(y2))/3;

plot(x1, y1, 'R', x2, y2)

grid on

%Получим следующие графики:

 

 

Листинг 1. Файл F.m.

function z=F(x,y)

z=sin(y+1)-x-1.2;

 

Листинг 2. Файл G.m.

function z=G(x,y)

z=2*y-cos(x)-3;

Листинг 3. Файл Fy.m.

function z=Fy(x,y)

z=cos(y+1);

 

Листинг 4. Файл Fx.m.

function z=Fx(x,y)

z=-1;

 

Листинг 5. Файл Gх_6.m.

function z=Gx(x,y)

z=sin(x);

Листинг 6. Файл Gy_6.m.

function z=Gy(x,y)

z=2;

Создайте файл Nuton.m (листинг 7), содержащий описание функции,

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

function Nuton(f,g,fx,fy,gx,gy,x0,y0,eps)

x1 = x0 +(feval(g,x0,y0)*feval(fy,x0,y0)-

feval(f,x0,y0)*feval(gy,x0,y0))/(feval(fx,x0,y0)*feval(gy,x0,y0)-feval(fy,x0,y0)*feval(gx,x0,y0));

y1=y0+(feval(f,x0,y0)*feval(gx,x0,y0)-feval(g,x0,y0)*feval(fx,x0,y0))/(feval(fx,x0,y0)*feval(gy,x0,y0)-feval(fy,x0,y0)*feval(gx,x0,y0));

k = 1;

while abs(x1-x0)>eps & abs(y1-y0)>eps

x0=x1;

y0=y1;

x1=x0+(feval(g,x0,y0)*feval(fy,x0,y0)-

feval(f,x0,y0)*feval(gy,x0,y0))/(feval(fx,x0,y0)*feval(gy,x0,y0)-feval(fy,x0,y0)*feval(gx,x0,y0));

y1=y0+(feval(f,x0,y0)*feval(gx,x0,y0)-feval(g,x0,y0)*feval(fx,x0,y0))/(feval(fx,x0,y0)*feval(gy,

x0,y0)-feval(fy,x0,y0)*feval(gx,x0,y0));

k=k+1;

end;

x=x1

y=y1

k

 

Найдите решение системы:

 

Nuton('F','G','Fx','Fy','Gx','Gy',-1,2,0.000001)

 

 

 

 

 

 

 

 

 

 

#include <iostream> 

 

using namespace std; 

 

#define eps 0.00001 

 

double function1(double x, double y)

{ 

   return sin(y + y) - 1.122*x - 0.498;

} 

 

double function2(double x, double y)

{ 

   return x*x + y*y - 1;

} 

 

double func11(double x, double y)

{ 

   return cos(x + y) - 1.122;

} 

 

double func12(double x, double y)

{ 

   return cos(x + y);

} 

 

double func21(double x, double y)

{ 

   return 2*x;

} 

 

double func22(double x, double y)

{ 

   return 2*y;

} 

 

void obr_matr(double a[2][2])

{ 

   double det, aa; 

   det = a[0][0]*a[1][1] - a[0][1]*a[1][0]; 

   aa = a[0][0]; 

   a[0][0] = a[1][1]/det; 

   a[1][1] = aa/det; 

   aa = a[0][1]; 

   a[0][1] = -a[1][0]/det; 

   a[1][0] = -aa/det;

} 

 

void nuton(double x, double y)

{ 

   int i = 1; 

   double a[2][2], dx, dy, b[2], norm; 

   do 

   { 

       a[0][0] = func11(x, y); 

       a[0][1] = func12(x, y); 

       a[1][0] = func21(x, y); 

       a[1][1] = func22(x, y); 

       ober_matr(a); 

       dx = -a[0][0]*function1(x, y) + -a[0][1]*function2(x, y); 

       dy = -a[1][0]*function1(x, y) + -a[1][1]*function2(x, y); 

       x = x + dx; 

       y = y + dy; 

       b[0] = function1(x, y); 

       b[1] = function2(x, y); 

       norm = sqrt(b[0]*b[0]+b[1]*b[1]); 

       i++; 

   } 

   while (norm >= eps);

cout << x << endl << y << endl;

}

void main()

{ 

   double x, y; 

   cout << "x = "; 

   cin >> x ; 

   cout << "y = "; 

   cin >> y; 

   nuton(x, y); 

   cout << endl; 

   system("PAUSE");

}

 

 

 

 

 

 

 

 

 

 

 

Задание 2. Решить систему уравнений упрощённым методом НЬЮТОНА

 с точностью ε=0,000001:

Система:

Перепишем систему в виде:

X=sin(y+1)-1,2

Y=(3+cos(x))/2

Графически находим  начальное приближение    x⁽⁰⁾=-1; y⁽⁰⁾=2

 

 

 

 

Перепишем систему в следующем  виде:

f₁(x,y)= sin(y+1)-x -1,2=0

f₂(x,y)=2y-cos x-3=0

Найдём частные производные:

(f₁)^’ _x (x,y)= -1;

(f₂ )^’_x(x,y)=sin(x);

-------------------------------------

(f₁)^’ _y(x,y)= cos(y+1)-;

(f₂ )^’_y(x,y)=2;

Матрица Якоби:

-1         cos(y+1)

sin(x)       2

 

Имеем СЛАУ:

-1∆x^(k)+cos(y+1)∆y^(k)= -f₁(x^(k),y^(k)),

sin(x^(k))∆x^(k)+2∆y^(k)= -f₂(x^(k),y^(k)

Коэффициенты  вычисляем 1 раз на 1 итерации :

 

-1∆x^(k)-0.9900∆y^(k)= -sin(y^(k)+1)+x^(k) +1.2,

-0.8415*∆x^(k)+2∆y^(k)= -2y^(k)+cos x^(k)+3

 

1 шаг полностью  совпадает с классическим методом

 

∆x⁽⁰⁾=0.1191

∆y⁽⁰⁾=-0,1798

x⁽¹⁾=x⁽⁰⁾+ ∆x⁽⁰⁾=-1+0.1191=-0,8809;

y⁽¹⁾=1,8202

2 шаг:

 Δy⁽¹⁾= -0,0024

    Δx⁽¹⁾= -0,0009

 

x⁽²⁾=x⁽¹⁾+ ∆x⁽¹⁾= -0,8818;

y⁽²⁾= 1,8179

3 шаг: 

Δy⁽²⁾= -0,000008

Δx⁽²⁾= -0,000091

x⁽³⁾=x⁽²⁾+ ∆x⁽²⁾= -0,881880;

y⁽³⁾= 1,817847

 

4 шаг:

Δy⁽³⁾= 0,00000214

Δx⁽³⁾= 0,00000246

x⁽⁴⁾ =-0,881883

y⁽⁴⁾ =1,817849

 

5 шаг:

Δy⁽⁴⁾==0,00000009

Δx⁽⁴⁾=0,000000003

 

x(5) =-0,881883

y(5) =1,8178495

 

Неравенство выполняется, поэтому итерации прекращаем.

ОТВЕТ: x =-0,8819

                  y =1,8178

 

 

MATLAB:

 

x1=-2:0.1:2;

y1=sin(x1-0.6)-1.6;

y2=-3:0.1:3;

x2=(0.9+cos(y2))/3;

plot(x1, y1, 'R', x2, y2)

grid on

%Получим следующие графики:

 

 

Листинг 1. Файл F.m.

function z=F(x,y)

z=sin(y+1)-x-1.2;

 

Листинг 2. Файл G.m.

function z=G(x,y)

z=2*y-cos(x)-3;

Листинг 3. Файл Fy.m.

function z=Fy(x,y)

z=cos(y+1);

 

Листинг 4. Файл Fx.m.

function z=Fx(x,y)

z=-1;

 

Листинг 5. Файл Gх_6.m.

function z=Gx(x,y)

z=sin(x);

Листинг 6. Файл Gy_6.m.

function z=Gy(x,y)

z=2;

Создайте файл Nuton.m (листинг 7), содержащий описание функции,

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

function Nuton(f,g,fx,fy,gx,gy,x0,y0,eps)

x1 = x0 +(feval(g,x0,y0)*feval(fy,x0,y0)-

feval(f,x0,y0)*feval(gy,x0,y0))/(feval(fx,x0,y0)*feval(gy,x0,y0)-feval(fy,x0,y0)*feval(gx,x0,y0));

y1=y0+(feval(f,x0,y0)*feval(gx,x0,y0)-feval(g,x0,y0)*feval(fx,x0,y0))/(feval(fx,x0,y0)*feval(gy,x0,y0)-feval(fy,x0,y0)*feval(gx,x0,y0));

k = 1;

while abs(x1-x0)>eps & abs(y1-y0)>eps

x0=x1;

y0=y1;

x1=x0+(feval(g,x0,y0)*feval(fy,x0,y0)-

feval(f,x0,y0)*feval(gy,x0,y0))/(feval(fx,x0,y0)*feval(gy,x0,y0)-feval(fy,x0,y0)*feval(gx,x0,y0));

y1=y0+(feval(f,x0,y0)*feval(gx,x0,y0)-feval(g,x0,y0)*feval(fx,x0,y0))/(feval(fx,x0,y0)*feval(gy,

x0,y0)-feval(fy,x0,y0)*feval(gx,x0,y0));

k=k+1;

end;

x=x1

y=y1

k

 

Найдите решение  системы:

 

Nuton('F','G','Fx','Fy','Gx','Gy',-1,2,0.000001)

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Задание 3 . сравнение 2 методов 

• Если метод Ньютона сходится, то сходимость имеет 2 порядок
• Условия сходимости:
• Начальное приближение достаточно близко к корню
• Функции f_i(x₁,x₂,…x_n) и их первые и вторые производные удовлетворяли некоторым ограничениям
• Недостаток:
• На каждом шаге надо вычислять n² частный производных
• (в методе простой итерации – всего n)
• Выход: модификация (упрощенный метод Ньютона)
• Фиксировать Якобиан, пересчитывать его лишь иногда

 

 

1) условия сходимости  метода Ньютона зависят

а) от значений первых и вторых производных функций невязок  по

искомым параметрам;

б) от близости предыдущего, а в конечном итоге начального

приближения к решению. При  этом погрешность последующего приближения

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

этом смысле говорят о квадратичной сходимости метода Ньютона.

2) по мере приближения  к решению сходимость резко  ускоряется.

3) при задании начального  приближения достаточно далеко  от решения

итерационный процесс  метода Ньютона может быть расходящимся.