MLoop 2

02/11/2015 by Женя Максимова

Задача. Используйте метод хорд для того, чтобы отыскать с точностью [latex]\varepsilon[/latex] все действительные корни уравнения [latex]\frac{x}{2 \cdot \sin x +1}=\tan(\ln(x^2+1))[/latex]. Для подготовки необходимых графиков воспользуйтесь этим ресурсом.

Тесты(найдено с помощью математической системы WolframAlpha):

[latex]A[/latex]

[latex]B[/latex]

[latex]x\approx[/latex]

-20

-11.6945378230838209122818536587051434153…

-1.25741503276862309237205903178504130394…

0.547316310185252929580383582338832450320…

10.9948442206261587135425985750810372810…

Код программы

#include <stdio.h>
#include <math.h>
#include <iostream>
#include<cfloat>
using namespace std;

const int max_iterations=100;
const double interval = 0.001;  // interval length for roots subdivision
bool is_root = false;  // is false if function has a second type break in root candidate point
double delta; //precision of root search
double A, B;

double f(double x)//our function
{
 double r = DBL_MAX;
 try
  {
   r =  x / (2*sin(x)+1) - tan(log(x*x+1));  
  }
 catch(...)
  {
   //nothing to do for arithmetic exception
  }
 return r;
}

bool is_change_sign(double fa, double fb)
{
 return fa*fb<0;
}

double find_root(double a, double b, double delta,  bool *yes)
{
 double fa;
 double fb;
 double fc;
 double c;
 int l=0;
 fa=f(a);
 fb=f(b);
 double dy,dy1; //old and new OY projections of chord
 dy=fabs(fa-fb);
 dy1=DBL_MAX;
 while(true)
   {
    l++;
    if(l > max_iterations)
      {
       *yes=false;
       return (a+b)/2;
      }
    c = a - (b - a) * fa/(fb - fa);
    fc = f(c);
    if(fa*fc>0)
      {
       a=c;
       fa=fc;
      }
    else
      {
       b=c;
       fb=fc;
      }
    if(fabs(fa)<delta)
      {
       *yes=true;
       return a;
      }
    if(fabs(fb)<delta)
      {
       *yes=true;
       return b;
      }
    if(l>10) //there's a magic number, we must check conditions of break or stop after some iterations(some = 10)
      {
       if(dy1>dy) //checking if there's break of second kind(pole)
         {
          *yes=false;
          return (a+b)/2;
         }
       if(fabs(b-a)<delta) //else we are near the root
         {
          *yes=true;
          return (a+b)/2;
         }
      }
    dy=dy1; //change old chord height
    dy1=fabs(fa-fb); //calculate new chord height
   }
 // we can't get here
 *yes=false;
 return c;
}

int main()
{
 scanf("%lf %lf %lf",&A,&B,&delta);
 double a,b;
 double fa, fb;
 /*
 test interval for finding roots
 int k=2;
 A= sqrt(exp(M_PI*(1/2.0+k))-1)-M_PI;//70010637.8772-10;
 B= sqrt(exp(M_PI*(1/2.0+k))-1)+M_PI;//70010637.8772+10;
 delta = 1e-12;
 printf("A=%20.19g A=%20.19g delta=%20.19g \n",A,B,delta);
 */
 fb = f(A);
 for(a = A, b = A + interval; b<=B; a += interval, b += interval)
  {
   fa = fb;
   fb = f(b);
   if(is_change_sign(fa,fb))
     {
      double r=find_root(a,b,delta, &is_root);
      if(is_root)
       {
        printf("%14.15lf \n",r);
       }
     }
  }
 return 0;
}

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

#include <stdio.h>

#include <math.h>

#include <iostream>

#include<cfloat>

using namespace std;

const int max_iterations=100;

const double interval = 0.001; // interval length for roots subdivision

bool is_root = false; // is false if function has a second type break in root candidate point

double delta; //precision of root search

double A, B;

double f(double x)//our function

{

double r = DBL_MAX;

try

{

r = x / (2*sin(x)+1) - tan(log(x*x+1));

}

catch(...)

{

//nothing to do for arithmetic exception

}

return r;

}

bool is_change_sign(double fa, double fb)

{

return fa*fb<0;

}

double find_root(double a, double b, double delta, bool *yes)

{

double fa;

double fb;

double fc;

double c;

int l=0;

fa=f(a);

fb=f(b);

double dy,dy1; //old and new OY projections of chord

dy=fabs(fa-fb);

dy1=DBL_MAX;

while(true)

{

l++;

if(l > max_iterations)

{

*yes=false;

return (a+b)/2;

}

c = a - (b - a) * fa/(fb - fa);

fc = f(c);

if(fa*fc>0)

{

a=c;

fa=fc;

}

else

{

b=c;

fb=fc;

}

if(fabs(fa)<delta)

{

*yes=true;

return a;

}

if(fabs(fb)<delta)

{

*yes=true;

return b;

}

if(l>10) //there's a magic number, we must check conditions of break or stop after some iterations(some = 10)

{

if(dy1>dy) //checking if there's break of second kind(pole)

{

*yes=false;

return (a+b)/2;

}

if(fabs(b-a)<delta) //else we are near the root

{

*yes=true;

return (a+b)/2;

}

dy=dy1; //change old chord height

dy1=fabs(fa-fb); //calculate new chord height

}

// we can't get here

*yes=false;

return c;

}

int main()

{

scanf("%lf %lf %lf",&A,&B,&delta);

double a,b;

double fa, fb;

test interval for finding roots

int k=2;

A= sqrt(exp(M_PI*(1/2.0+k))-1)-M_PI;//70010637.8772-10;

B= sqrt(exp(M_PI*(1/2.0+k))-1)+M_PI;//70010637.8772+10;

delta = 1e-12;

printf("A=%20.19g A=%20.19g delta=%20.19g \n",A,B,delta);

fb = f(A);

for(a = A, b = A + interval; b<=B; a += interval, b += interval)

{

fa = fb;

fb = f(b);

if(is_change_sign(fa,fb))

{

double r=find_root(a,b,delta, &is_root);

if(is_root)

{

printf("%14.15lf \n",r);

}

return 0;

}

Алгоритм

Для начала запишем данное нам уравнение в виде функции [latex]y=f(x)[/latex] и построим ее график:

[latex]y=\frac{x}{2 \cdot sin x +1}-tan(ln(x^2+1))[/latex]

Задача о нахождении приближённых значений действительных корней уравнения [latex]f(x)=0[/latex] предусматривает предварительное отделение корня, то есть установление интервала, в котором других корней данного уравнения нет.

Метод хорд предусматривает замену функции на интервале на секущую, и поиск ее пересечения с осью [latex]OX[/latex]. На заданном интервале [latex][a,b][/latex] с точностью [latex]\varepsilon[/latex] корень будет вычисляться согласно итерационному выражению, которое имеет вид:

[latex]x_{i+1}=x_{i-1}-\frac{f(x_{i-1}) \cdot (x_{i}-x_{i-1})}{f(x_{i})-f(x_{i-1}) ) }[/latex]

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

В точке, где находится предполагаемый корень, имеется разрыв второго рода. Здесь метод хорд обнаруживает перемену знака и начинает сужать отрезок. Однако расстояние между крайними точками отрезка не уменьшается, а увеличивается. А именно увеличивается проекция отрезка на ось [latex]OY[/latex]. И чем ближе мы находимся к точке разрыва, тем она больше, а в самой точке стремится к бесконечности. В качестве примера приведем функцию [latex]\frac{(x-5)^2}{x-4}[/latex], имеющую разрыв второго рода в точке [latex]x=4[/latex].
Аналогичный случай — точка разрыва первого рода, где наша хорда стремится к некоторой константе — величине разрыва. Пример — функция [latex]\frac{\sin x\cdot(x-2.5)}{\left | x-2,5 \right |}[/latex], имеющая разрыв первого рода в точке [latex]x=2,5[/latex].
Функция не является разрывной и даже имеет корень, однако в точке корня производная стремится к бесконечности. Например, функция [latex]\sqrt[3]{x-1,2}[/latex], имеющая корень в точке [latex]x=1,2[/latex]. Для функций подобного рода длина проекции отрезка на ось [latex]OY[/latex] будет очень медленно меняться, вследствие чего для разумного числа итераций она будет превосходить заранее выбранную точность [latex]\varepsilon[/latex].

Функция равна нулю или очень близка к нулю на некотором интервале(например, функция [latex]y=rect(x-1,5)\cdot(x-1)\cdot(x-2)^2[/latex]). Здесь метод хорд найдет пересечение с осью [latex]OX[/latex] в интервале, где находятся корни(одна из сторон отрезка будет корнем), но все последующие итерации будут выдавать эту же точку. Поэтому хорда не будет уменьшаться, и даже этот один корень не будет найден(если будет использоваться стандартная [latex]\delta[/latex]-оценка точности по оси [latex]OX[/latex]).

Функция вида [latex]y=x^{2k}[/latex] или ей подобная, например, [latex]y=1+\sin x[/latex], к которой метод хорд вообще не применим, так как нарушается начальное условие применимости этого метода. Здесь нужно ввести дополнительную проверку. Изменим значение функции на небольшую константу — нашу точность [latex]\varepsilon[/latex] и повторим процедуру поиска корней. В результате мы получим [latex]2k[/latex] корней, каждую пару из которых мы можем считать краями интервала, в котором лежит настоящий корень.

К счастью, в данной нам функции присутствует только один из этих случаев, а именно разрыв второго рода. Аналитически рассмотрев нашу функцию, мы обнаружили, что корни следует искать в окрестности точек[latex]\sqrt{e^{\frac{\pi}{2}+\pi \cdot k}-1}[/latex] с отклонением [latex]\pm\pi[/latex]. В корнях функции ее производная быстро растет с ростом [latex]k[/latex].

Критерием отбрасывания кандидата на корень будет рост длины хорды при сужении интервала. Критерием останова будет сужение интервала до заданной точности [latex]\delta[/latex].
Код программы

Related Images:

Ю11.16

24/10/2014 by Іванов Вячеслав Володимирович

Задача:
Для заданной матрицы [latex]A(n, n)[/latex] найти обратную [latex]A^{-1}[/latex], используя итерационную формулу:
[latex]A_{k}^{-1} = A_{k-1}^{-1}(2E-A A_{k}^{-1}),[/latex] где [latex]E[/latex] — единичная матрица, [latex]A_{0}^{-1} = E[/latex]. Итерационный процесс заканчивается, если для заданной погрешности [latex]\varepsilon[/latex] справедливо:
[latex]\left| det(A A_{k}^{-1})-1 \right| \le \varepsilon[/latex]

Анализ задачи:
Прежде чем приступать к решению средствами языка C++, я создал прототип в системе компьютерной алгебры Wolfram Mathematica, с которым планировал сверяться при тестировании программы. Тем не менее, внимательный анализ показал, что с таким выбором начального приближения процесс уже на пятом шаге в значительной мере расходится даже для матриц размера 2*2. После уточнения условий и анализа дополнительного материала, посвященного методу Ньютона-Шульца, исходное приближение было изменено (по результатам исследования «An Improved Newton Iteration for the Generalized Inverse of a Matrix, with Applications», Victor Pan, Robert Schreiber, стр. 8):
[latex]A_{0}^{-1} =\frac { A }{ \left\| A \right\|_{1} \left\| A \right\| _{\infty } }[/latex], где [latex]{ \left\| A \right\| }_{ 1 }=\max _{ i }{ \sum _{ j=0 }^{ n-1 }{ \left| { a }_{ ij } \right| } } [/latex], [latex]{ \left\| A \right\| }_{ \infty }=\max _{ j }{ \sum _{ i=0 }^{ n-1 }{ \left| { a }_{ ij } \right| } }[/latex].
Эффективность предложенного подхода иллюстрируют результаты работы прототипа:

Следовательно, из пространства задачи можно переместиться в пространство решения и составить алгоритм реализации предложенного метода на языке C++.

Тесты:

[latex]n[/latex]	[latex]A[/latex]	[latex]A^{-1}[/latex]	Результат
3	1 2 3 5 5 7 11 13 7	-0.964771 0.430661 -0.017183 0.723533 -0.447973 0.137884 0.172358 0.155200 -0.086211	Пройден
2	1 2 2 1	-0.333067 0.666400 0.666400 -0.333067	Пройден
5	1 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 1	1 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 1	Пройден
3	1 2 3 4 5 6 7 8 9	Матрица вырождена	Пройден

Программный код:

#include <cstdio>
#include <algorithm>
#include <cmath>
//очистить выделенную память
void clear(double **arr, int n)
{
    for(int i = 0; i < n; i++)
        delete[] arr[i];
    delete[] arr;
}
//создать копию массива
double** clone(double **arr, int n)
{
    double **newArr = new double*[n];
    for(int row = 0; row < n; row++)
    {
        newArr[row] = new double[n];
        for(int col = 0; col < n; col++)
            newArr[row][col] = arr[row][col];
    }
    return newArr;
}
//print the matrix
void show(double **matrix, int n)
{
    for(int row = 0; row < n; row++){
            for(int col = 0; col < n; col++)
                printf("%lf\t", matrix[row][col]);
            printf("\n");
    }
    printf("\n");
}
//матричное умножение матриц
double** matrix_multi(double **A, double **B, int n)
{
    double **result = new double*[n];
    //заполнение нулями
    for(int row = 0; row < n; row++){
        result[row] = new double[n];
        for(int col = 0; col < n; col++){
            result[row][col] = 0;
        }
    }

    for(int row = 0; row < n; row++){
        for(int col = 0; col < n; col++){
            for(int j = 0; j < n; j++){
                result[row][col] += A[row][j] * B[j][col];
            }
        }
    }
    return result;
}
//умножение матрицы на число
void scalar_multi(double **m, int n, double a){
    for(int row = 0; row < n; row++)
        for(int col = 0; col < n; col++){
            m[row][col] *= a;
        }
}
//вычисление суммы двух квадратных матриц
void sum(double **A, double **B, int n)
{
    for(int row = 0; row < n; row++)
        for(int col = 0; col < n; col++)
            A[row][col] += B[row][col];
}

//вычисление определителя
double det(double **matrix, int n) //квадратная матрица размера n*n
{
    double **B = clone(matrix, n);
    //приведение матрицы к верхнетреугольному виду
    for(int step = 0; step < n - 1; step++)
        for(int row = step + 1; row < n; row++)
        {
            double coeff = -B[row][step] / B[step][step]; //метод Гаусса
            for(int col = step; col < n; col++)
                B[row][col] += B[step][col] * coeff;
        }
    //Рассчитать определитель как произведение элементов главной диагонали
    double Det = 1;
    for(int i = 0; i < n; i++)
        Det *= B[i][i];
    //Очистить память
    clear(B, n);
    return Det;
}

int main()
{
    //Исходная матрица, динамический двухмерный массив
    int n; scanf("%d", &n);
    double **A = new double*[n];
    for(int row = 0; row < n; row++)
    {
        A[row] = new double[n];
            for(int col = 0; col < n; col++)
                scanf("%lf", &A[row][col]);
    }

    /* Численное вычисление обратной матрицы по методу Ньютона-Шульца
        1. Записать начальное приближение [Pan, Schreiber]:
            1) Транспонировать данную матрицу
            2) Нормировать по столбцам и строкам
        2. Повторять процесс до достижения заданной точности.
    */

    double N1 = 0, Ninf = 0; //норма матрицы по столбцам и по строкам
    double **A0 = clone(A, n);       //инициализация начального приближения
    for(size_t row = 0; row < n; row++){
        double colsum = 0, rowsum = 0;
        for(size_t col = 0; col < n; col++){
            rowsum += fabs(A0[row][col]);
            colsum += fabs(A0[col][row]);
        }
        N1 = std::max(colsum, N1);
        Ninf = std::max(rowsum, Ninf);
    }
    //транспонирование
    for(size_t row = 0; row < n - 1; row++){
        for(size_t col = row + 1; col < n; col++)
            std::swap(A0[col][row], A0[row][col]);
    }
    scalar_multi(A0, n, (1 / (N1 * Ninf))); //нормирование матрицы
    //инициализация удвоенной единичной матрицы нужного размера
    double **E2 = new double*[n];
    for(int row = 0; row < n; row++)
    {
        E2[row] = new double[n];
            for(int col = 0; col < n; col++){
                if(row == col)
                    E2[row][col] = 2;
                else
                    E2[row][col] = 0;
            }
    }
    double **inv = clone(A0, n); //A_{0}
    double EPS = 0.001;   //погрешность
    if(det(A, n) != 0){ //если матрица не вырождена
        while(fabs(det(matrix_multi(A, inv, n), n) - 1) >= EPS) //пока |det(A * A[k](^-1)) - 1| >= EPS
        {
            double **prev = clone(inv, n); //A[k-1]
            inv = matrix_multi(A, prev, n);   //A.(A[k-1]^(-1))
            scalar_multi(inv, n, -1);         //-A.(A[k-1]^(-1))
            sum(inv, E2, n);                   //2E - A.(A[k-1]^(-1))
            inv = matrix_multi(prev, inv, n); //(A[k-1]^(-1)).(2E - A.(A[k-1]^(-1)))
            clear(prev, n);
        }
        //вывод матрицы на экран
        show(inv, n);
    }
    else
        printf("Impossible\n");
    clear(A, n);
    clear(E2, n);
    return 0;
}

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

126

127

128

129

130

131

132

133

134

135

136

137

138

139

140

141

142

143

144

145

146

147

148

149

150

151

152

153

154

155

156

157

158

#include <cstdio>

#include <algorithm>

#include <cmath>

//очистить выделенную память

void clear(double **arr, int n)

{

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

delete[] arr[i];

delete[] arr;

}

//создать копию массива

double** clone(double **arr, int n)

{

double **newArr = new double*[n];

for(int row = 0; row < n; row++)

{

newArr[row] = new double[n];

for(int col = 0; col < n; col++)

newArr[row][col] = arr[row][col];

}

return newArr;

}

//print the matrix

void show(double **matrix, int n)

{

for(int row = 0; row < n; row++){

for(int col = 0; col < n; col++)

printf("%lf\t", matrix[row][col]);

printf("\n");

}

printf("\n");

}

//матричное умножение матриц

double** matrix_multi(double **A, double **B, int n)

{

double **result = new double*[n];

//заполнение нулями

for(int row = 0; row < n; row++){

result[row] = new double[n];

for(int col = 0; col < n; col++){

result[row][col] = 0;

}

for(int row = 0; row < n; row++){

for(int col = 0; col < n; col++){

for(int j = 0; j < n; j++){

result[row][col] += A[row][j] * B[j][col];

}

return result;

}

//умножение матрицы на число

void scalar_multi(double **m, int n, double a){

for(int row = 0; row < n; row++)

for(int col = 0; col < n; col++){

m[row][col] *= a;

}

//вычисление суммы двух квадратных матриц

void sum(double **A, double **B, int n)

{

for(int row = 0; row < n; row++)

for(int col = 0; col < n; col++)

A[row][col] += B[row][col];

}

//вычисление определителя

double det(double **matrix, int n) //квадратная матрица размера n*n

{

double **B = clone(matrix, n);

//приведение матрицы к верхнетреугольному виду

for(int step = 0; step < n - 1; step++)

for(int row = step + 1; row < n; row++)

{

double coeff = -B[row][step] / B[step][step]; //метод Гаусса

for(int col = step; col < n; col++)

B[row][col] += B[step][col] * coeff;

}

//Рассчитать определитель как произведение элементов главной диагонали

double Det = 1;

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

Det *= B[i][i];

//Очистить память

clear(B, n);

return Det;

}

int main()

{

//Исходная матрица, динамический двухмерный массив

int n; scanf("%d", &n);

double **A = new double*[n];

for(int row = 0; row < n; row++)

{

A[row] = new double[n];

for(int col = 0; col < n; col++)

scanf("%lf", &A[row][col]);

}

/* Численное вычисление обратной матрицы по методу Ньютона-Шульца

1. Записать начальное приближение [Pan, Schreiber]:

1) Транспонировать данную матрицу

2) Нормировать по столбцам и строкам

2. Повторять процесс до достижения заданной точности.

double N1 = 0, Ninf = 0; //норма матрицы по столбцам и по строкам

double **A0 = clone(A, n); //инициализация начального приближения

for(size_t row = 0; row < n; row++){

double colsum = 0, rowsum = 0;

for(size_t col = 0; col < n; col++){

rowsum += fabs(A0[row][col]);

colsum += fabs(A0[col][row]);

}

N1 = std::max(colsum, N1);

Ninf = std::max(rowsum, Ninf);

}

//транспонирование

for(size_t row = 0; row < n - 1; row++){

for(size_t col = row + 1; col < n; col++)

std::swap(A0[col][row], A0[row][col]);

}

scalar_multi(A0, n, (1 / (N1 * Ninf))); //нормирование матрицы

//инициализация удвоенной единичной матрицы нужного размера

double **E2 = new double*[n];

for(int row = 0; row < n; row++)

{

E2[row] = new double[n];

for(int col = 0; col < n; col++){

if(row == col)

E2[row][col] = 2;

else

E2[row][col] = 0;

}

double **inv = clone(A0, n); //A_{0}

double EPS = 0.001; //погрешность

if(det(A, n) != 0){ //если матрица не вырождена

while(fabs(det(matrix_multi(A, inv, n), n) - 1) >= EPS) //пока |det(A * A[k](^-1)) - 1| >= EPS

{

double **prev = clone(inv, n); //A[k-1]

inv = matrix_multi(A, prev, n); //A.(A[k-1]^(-1))

scalar_multi(inv, n, -1); //-A.(A[k-1]^(-1))

sum(inv, E2, n); //2E - A.(A[k-1]^(-1))

inv = matrix_multi(prev, inv, n); //(A[k-1]^(-1)).(2E - A.(A[k-1]^(-1)))

clear(prev, n);

}

//вывод матрицы на экран

show(inv, n);

}

else

printf("Impossible\n");

clear(A, n);

clear(E2, n);

return 0;

}

Программа доступна для тестирования по ссылке: http://ideone.com/7YphoX.

Алгоритм решения
При решении данной задачи оправдывает себя подход «разделяй и властвуй»: вычисление обратной матрицы подразумевает промежуточные этапы, корректной реализации которых следует уделить особое внимание. Наверняка стоило бы написать класс matrix и реализовать перегрузку операций, но задачу можно решить и без применения объектно-ориентированных средств, пусть от этого решение и потеряет в изящности.
Можно выделить следующие подзадачи:

Инициализация динамического двухмерного массива и передача его в функцию.
Вычисление определителя матрицы (с применением метода Гаусса).
Вычисление нормы матрицы.
Транспонирование матрицы.
Умножение матрицы на скаляр.
Матричное умножение матриц.
Сложение двух матриц.
Непосредственно итерационный процесс Ньютона-Шульца

Ниже приведено пояснение к подходу к реализации некоторых подзадач:

Выделение памяти для хранения массива происходит не на этапе запуска программы, после компиляции. Для инициализации использован конструктор new
Вычисление определителя можно разбить на два последовательных шага:
1. Приведение матрицы к верхнетреугольному виду (методом Гаусса).
2. Вычисление определителя как произведения элементов главной диагонали.
Нормы матрицы вычисляются как максимальные значения суммы элементов по столбцам и строкам соответственно.
При транспонировании обмениваются местами элементы, симметричные главной диагонали.
При умножении матрицы на скаляр каждый элемент матрицы умножается на соответствующее вещественное число.
При перемножении двух квадратных матриц используется промежуточный массив для хранения результата вычислений.
Сложение двух матриц аналогично попарному сложению элементов, расположенных на соответствующих позициях.
Максимально допустимая погрешность для метода Ньютона-Шульца [latex]\varepsilon = 0.001[/latex]. Программа использует локальный массив для хранения [latex]A_{k}^{-1}[/latex], инициализация которого происходит в теле цикла.

Технические детали реализации:
При выполнении подзадач часто приходится использовать локальные массивы, так что для очистки выделенной под них памяти создана отдельная функция clear().