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

№10

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

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

x₂

0.5 x₁

Метод Ньютона решения системы состоит в построении итерационной последовательности:

Запишем в векторном виде.

,

Из графика

Составим матрицу Якоби.

- 0.9098 – 0.3768 = -1.2866 можно составить обратную матрицу.

,

По формуле (1) получаем

Т.о.

Найдем .

Т.о.

Алгоритм на C++

#include <iostream>

#include <cmath>

#include <vector>

using namespace std;

double det(vector <vector <double> > a)

{

return a[0][0]*a[1][1] - a[0][1]*a[1][0];

}

vector <vector <double>> inverse(vector <vector <double>> a)

{

vector <vector <double>> inv(a.size(), vector <double> (a.size()));

double d = det(a);

inv[0][1] = -a[0][1];

inv[1][0] = -a[1][0];

inv[0][0] = a[1][1];

inv[1][1] = a[0][0];

for (int i = 0; i < inv.size(); i++)

for (int j = 0; j < inv.size(); j++)

inv[i][j] /= d;

return inv;

}

vector <vector <double>> Jacobian(double x1, double x2)

{

vector <vector <double>> J(2, vector <double> (2));

J[0][0] = 2 * x1;

J[0][1] = 2 * x2;

J[1][0] = 2 * x1;

J[1][1] = -1;

return J;

}

double f1(double x1, double x2)

{

return x1*x1 + x2*x2 - 0.25;

}

double f2(double x1, double x2)

{

return x1*x1 - x2;

}

vector <double> f(double x1, double x2)

{

vector <double> g(2);

g[0] = f1(x1, x2);

g[1] = f2(x1, x2);

return g;

}

vector <double> x (double x1, double x2)

{

vector <double> g(2);

g[0] = x1;

g[1] = x2;

return g;

}

vector <double> operator * (vector <vector <double>> a, vector <double> b)

{

vector <double> res(2);

for (int i = 0; i < 2; i++)

for (int j = 0; j < 2; j++)

res[i] += a[i][j] * b[j];

return res;

}

vector <double> operator - (vector <double> a, vector <double> b)

{

vector <double> res(2);

for (int i = 0; i < 2; i++)

res[i] = a[i] - b[i];

return res;

}

int main()

{

freopen ("input.txt", "rt", stdin);

freopen ("output.txt", "wt", stdout);

double EPS = 1E-9;

double x1 = 0.4549;

double x2 = 0.20709;

while (true)

{

vector <double> xn = x(x1, x2) - inverse(Jacobian(x1, x2))*f(x1, x2);

if (fabs (xn[0] - x1) < EPS && fabs (xn[1] - x2) < EPS) break;

x1 = xn[0];

x2 = xn[1];

}

printf("%.11lf %.11lf", x1, x2);

return 0;

}

Дата добавления: 2015-08-28