Метод Ньютона - Общие вопросы C/C++

Symple me · 08.11.2012, 12:25

Найдите, пожалуйста, ошибку!

Код:

#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <iostream>
#include <conio.h>
#pragma hdrstop
using namespace std;
//---------------------------------------------------------------------------
float max(float m[4])
{
float max=fabs(m[0]);
for (int i=0; i<4; i++)
if (fabs(m[i])>max)
max=fabs(m[i]) ;
return max;
}
void inversion(float A[4][4], int N)         //функция, находящая
{                                            //обратную матрицу
    float temp;
    float B[4][4];
    float res[4][4];
    int i,j,k;

    float **E = new float *[N];

    for (int i = 0; i < N; i++)
        E[i] = new float [N];
        for (i=0; i<N; i++)
        for (j=0; j<N; j++)
        B[i][j]=A[i][j];

    for (int i = 0; i < N; i++)
        for (int j = 0; j < N; j++)
        {
            E[i][j] = 0.0;

            if (i == j)
                E[i][j] = 1.0;
        }

    for (int k = 0; k < N; k++)
    {
        temp = A[k][k];

        for (int j = 0; j < N; j++)
        {
            A[k][j] /= temp;
            E[k][j] /= temp;
        }

        for (int i = k + 1; i < N; i++)
        {
            temp = A[i][k];

            for (int j = 0; j < N; j++)
            {
                A[i][j] -= A[k][j] * temp;
                E[i][j] -= E[k][j] * temp;
            }
        }
    }

    for (int k = N - 1; k > 0; k--)
    {
        for (int i = k - 1; i >= 0; i--)
        {
            temp = A[i][k];

            for (int j = 0; j < N; j++)
            {
                A[i][j] -= A[k][j] * temp;
                E[i][j] -= E[k][j] * temp;
            }
        }
    }

    for (int i = 0; i < N; i++)
        for (int j = 0; j < N; j++)
            A[i][j] = E[i][j];
                 for (int i = 0; i < N; i++)
        delete [] E[i];

    delete [] E;
}
float J00( float a, float b)                 //функции, возвращающие
{                                            //производные 4 уравнений
float J00 =(35*b*sin(b+7*a))/16;             //системы по 4 неизвестным
return J00;
}
float J01(float a, float b)
{
float J01=-(cos(b+7*a)+b*sin(b+7*a))/32;
return J01;
}
float J02(float c)
{
float J02=cos(c-3)/16;
return J02;
}
float J03()
{
return 0;
}
float J10 (float a, float b, float c)
{
float J10=(c+4-(b/(2*sqrt(a))))/256;
return J10;
}
float J11(float a, float c)
{
float J11=-((3*c+12+sqrt(a))/256);
return J11;
}
float J12 (float a, float b)
{
float J12=(a-3*b)/256;
return J12;
}
float J13 ()
{
return 0;
}
float J20(float b)
{
float J20=b/4096;
return J20;
}
float J21 (float a, float b, float c)
{
float J21=(a/4096)- c*cos(c*b + 0.5);
return J21;
}
float J22 (float b, float c)
{
float J22=-cos (c*b+0.5)*b;
return J22;
}
float J23 ()
{
return 0;
}
float J30 (float a, float d)
{
float J30=-d*sin(a*d-1.5)+(1/256);
return J30;
}
float J31 ()
{
return 0;
}
float J32 ()
{
return 0;
}
float J33 (float a, float d)
{
float J33=-a*sin (a*d-1.5);
return J33;
}
int main
()
{
const int N=4;
int i,j,k;
float J[N][N], Jinv[N][N];
float x[N],b[N],f[N];
float E, d;
E = 0.0000001;
k= 4;
printf("Vvedite priblizjenie:\n");
for (i=0; i<k; i++)
{
printf("x[%d]= ",i);
scanf("%f",&x[i]);
}
do
{
J[0][0]=J00(x[0],x[1]);                 //Матрица Якоби
J[0][1]=J01(x[0],x[1]);
J[0][2]=J02(x[2]);
J[0][3]=J03();
J[1][0]=J10 (x[0], x[1], x[2]);
J[1][1]=J11 (x[0], x[2]);
J[1][2]=J12 (x[0], x[1]);
J[1][3]=J13 ();
J[2][0]=J20 (x[1]);
J[2][1]=J21 (x[0], x[1], x[2]);
J[2][2]=J22 (x[1], x[2]);
J[2][3]=J23 ();
J[3][0]=J30 (x[0], x[3]);
J[3][1]=J31 ();
J[3][2]=J32 ();
J[3][3]=J33 (x[0], x[3]);
for (i=0; i<N; i++)
for (j=0; j<N; j++)
Jinv[i][j]=J[i][j];
inversion (Jinv, N);                          //Обратная матрица Якоби
f[0]=sin (x[2]-3)-0.5*x[1]*cos(x[1]+7*x[0])/16;     //значение функций
f[1]=((x[0]-3*x[1])*(x[2]+4)-sqrt(x[0])*x[1])/256; //от приближения
f[2]=((x[0]*x[1])/4096)-sin (x[2]*x[1]+0.5);
f[3]=cos (x[3]*x[0]-1.5)+x[0]/256;
for (int j=0;j<N;j++)
{
b[j]=0;
for (int i=0; i<N; i++)
b[j]=b[j]+Jinv[i][j]*f[i];                        //вектор поправок
}
for (i=0; i<N; i++)
{
x[i]=x[i]+b[i];
}
d=max(b);
}
while (d>=E);
for (i=0; i<k; i++)
printf("x[%d]= %8.3f",i, x[i]);
getch();
}

Murashov · 08.11.2012, 12:35

ужас.. там прога на строк 20-30 вроде максимум

Symple me · 08.11.2012, 16:27

Просто я вручную посчитала производные для каждого из 4 уравнений по каждой из 4 переменных. И с 86 по 160 строку я описываю 16 функций, возвращающие эти производные. Еще функция для нахождения обратной матрицы длинная: с 18 по 85 строку.

Murashov · 08.11.2012, 18:18

какие данные вводите?

Symple me · 08.11.2012, 19:56

Разные пробовала. Нули только не ввожу, потому что в некоторых производных встречается деление на компоненты начального приближения. Например, вводила 0.5 для каждого икса. Пишет в консоли:
sqrt: DOMAIN error.
А Билдер пишет: Invalid floating point operation.

Symple me · 13.11.2012, 06:14

Запустите, пожалуйста!

Код:

#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <iostream>
#include <conio.h>
#pragma hdrstop
using namespace std;
//---------------------------------------------------------------------------
void inversion(float A[4][4], int N)         //функция, для нахождения
{                                                    //обратной матрицы  
    float temp;
    float B[4][4];
    float res[4][4];
    int i,j,k;

    float **E = new float *[N];

    for (int i = 0; i < N; i++)
        E[i] = new float [N];
        for (i=0; i<N; i++)
        for (j=0; j<N; j++)
        B[i][j]=A[i][j];

    for (int i = 0; i < N; i++)
        for (int j = 0; j < N; j++)
        {
            E[i][j] = 0.0;

            if (i == j)
                E[i][j] = 1.0;
        }

    for (int k = 0; k < N; k++)
    {
        temp = A[k][k];

        for (int j = 0; j < N; j++)
        {
            A[k][j] /= temp;
            E[k][j] /= temp;
        }

        for (int i = k + 1; i < N; i++)
        {
            temp = A[i][k];

            for (int j = 0; j < N; j++)
            {
                A[i][j] -= A[k][j] * temp;
                E[i][j] -= E[k][j] * temp;
            }
        }
    }

    for (int k = N - 1; k > 0; k--)
    {
        for (int i = k - 1; i >= 0; i--)
        {
            temp = A[i][k];

            for (int j = 0; j < N; j++)
            {
                A[i][j] -= A[k][j] * temp;
                E[i][j] -= E[k][j] * temp;
            }
        }
    }

    for (int i = 0; i < N; i++)
        for (int j = 0; j < N; j++)
            A[i][j] = E[i][j];
                 for (int i = 0; i < N; i++)
        delete [] E[i];

    delete [] E;
}
float J00( float a, float b)                   //функции, возвращающие
{                                                   //производные 4 уравнений системы
float J00 =(35*b*sin(b+7*a))/16;       //по 4 неизвестным
return J00;
}
float J01(float a, float b)
{
float J01=-(cos(b+7*a)+b*sin(b+7*a))/32;
return J01;
}
float J02(float c)
{
float J02=cos(c-3)/16;
return J02;
}
float J03()
{
return 0;
}
float J10 (float a, float b, float c)
{
float J10=(c+4-(b/(2*sqrt(a))))/256;
return J10;
}
float J11(float a, float c)
{
float J11=-((3*c+12+sqrt(a))/256);
return J11;
}
float J12 (float a, float b)
{
float J12=(a-3*b)/256;
return J12;
}
float J13 ()
{
return 0;
}
float J20(float b)
{
float J20=b/4096;
return J20;
}
float J21 (float a, float b, float c)
{
float J21=(a/4096)- c*cos(c*b + 0.5);
return J21;
}
float J22 (float b, float c)
{
float J22=-cos (c*b+0.5)*b;
return J22;
}
float J23 ()
{
return 0;
}
float J30 (float a, float d)
{
float J30=-d*sin(a*d-1.5)+(1/256);
return J30;
}
float J31 ()
{
return 0;
}
float J32 ()
{
return 0;
}
float J33 (float a, float d)
{
float J33=-a*sin (a*d-1.5);
return J33;
}
int main()
{
const int N=4;
int i,j,k;
float J[N][N], Jinv[N][N];
float x[N],b[N],f[N];
float E;
E = 0.0000001;
printf("Vvedite priblizjenie:\n");
for (i=0; i<4; i++)
{
printf("x[%d]= ",i);
scanf("%f",&x[i]);
}
do
{
J[0][0]=J00(x[0],x[1]);                 //матрица Якоби
J[0][1]=J01(x[0],x[1]);
J[0][2]=J02(x[2]);
J[0][3]=J03();
J[1][0]=J10 (x[0], x[1], x[2]);
J[1][1]=J11 (x[0], x[2]);
J[1][2]=J12 (x[0], x[1]);
J[1][3]=J13 ();
J[2][0]=J20 (x[1]);
J[2][1]=J21 (x[0], x[1], x[2]);
J[2][2]=J22 (x[1], x[2]);
J[2][3]=J23 ();
J[3][0]=J30 (x[0], x[3]);
J[3][1]=J31 ();
J[3][2]=J32 ();
J[3][3]=J33 (x[0], x[3]);
for (i=0; i<N; i++)
for (j=0; j<N; j++)
Jinv[i][j]=J[i][j];
inversion (Jinv, N);                          //обратная матрица Якоби
f[0]=sin (x[2]-3)-0.5*x[1]*cos(x[1]+7*x[0])/16;     //значение функций
f[1]=((x[0]-3*x[1])*(x[2]+4)-sqrt(x[0])*x[1])/256; //от приближения
f[2]=((x[0]*x[1])/4096)-sin (x[2]*x[1]+0.5);
f[3]=cos (x[3]*x[0]-1.5)+x[0]/256;
for (int j=0;j<N;j++)
{
b[j]=0;
for (int i=0; i<N; i++)
b[j]=b[j]+Jinv[i][j]*f[i];                        //вектор поправок
}
for (i=0; i<N; i++)
x[i]=x[i]-b[i];                                      //уточнение решения
}
while (fabs(f[i])>E);
for (i=0; i<4; i++)
printf("x[%d]= %8.3f",i, x[i]);
getch();
}

Symple me · 13.11.2012, 06:21

Задается начальное приближение x[i]. На каждой итерации считается вектор поправок b[i] (в том числе на каждой итерации пересчитываются производные). Решение уточняется: x[i]=x[i]-b[i]. Некоторые компоненты x становятся отрицательными, а при подсчете производных из них извлекается квадратный корень. Поэтому числа косячные(

08.11.2012, 16:27	#3
Symple me Пользователь Регистрация: 09.10.2012 Сообщений: 47	Просто я вручную посчитала производные для каждого из 4 уравнений по каждой из 4 переменных. И с 86 по 160 строку я описываю 16 функций, возвращающие эти производные. Еще функция для нахождения обратной матрицы длинная: с 18 по 85 строку. There are 10 types of people: those who understand binary and those who don't.

08.11.2012, 19:56	#5
Symple me Пользователь Регистрация: 09.10.2012 Сообщений: 47	Разные пробовала. Нули только не ввожу, потому что в некоторых производных встречается деление на компоненты начального приближения. Например, вводила 0.5 для каждого икса. Пишет в консоли: sqrt: DOMAIN error. А Билдер пишет: Invalid floating point operation. There are 10 types of people: those who understand binary and those who don't.

13.11.2012, 06:21	#7
Symple me Пользователь Регистрация: 09.10.2012 Сообщений: 47	Задается начальное приближение x[i]. На каждой итерации считается вектор поправок b[i] (в том числе на каждой итерации пересчитываются производные). Решение уточняется: x[i]=x[i]-b[i]. Некоторые компоненты x становятся отрицательными, а при подсчете производных из них извлекается квадратный корень. Поэтому числа косячные( There are 10 types of people: those who understand binary and those who don't.

Похожие темы
Тема	Автор	Раздел	Ответов	Последнее сообщение
Метод Ньютона	S1LenseR	Помощь студентам	0	30.04.2012 15:09
Метод Ньютона	fariou73	Помощь студентам	1	24.04.2012 16:45
Методы оптимизации: метод Ньютона и метод наискорейшего спуска	ruslanGacurap	Помощь студентам	0	30.01.2012 13:54
Метод Ньютона	lordsyrius	Помощь студентам	5	23.11.2009 23:52
Метод Ньютона	Cubar	Помощь студентам	12	09.02.2008 21:28

08.11.2012, 12:35	#2
Murashov Форумчанин Регистрация: 30.10.2012 Сообщений: 121	ужас.. там прога на строк 20-30 вроде максимум

08.11.2012, 18:18	#4
Murashov Форумчанин Регистрация: 30.10.2012 Сообщений: 121	какие данные вводите?