Метод Ньютона в решении СНАУ

#include <math.h>
#include <conio.h>
#include <stdio.h>
#include <locale.h>
 
// Количество уравнений
const int n = 2;
 
void GaussSolve (float a[n][n],float b[n],float x[n]);
float f1(float x[n]);
float f2(float x[n]);
float df1_dx(float x[n]);
float df1_dy(float x[n]);
float df2_dx(float x[n]);
float df2_dy(float x[n]);
 
float norm (float x[n]);
 
int main()
{
setlocale (LC_ALL,"");
printf ("%s","Решение системы нелинейных уравнений методом Ньютона\n");
 
int i, iter=0;
 
// Точность
float eps = 0.00001f;
 
// Матрица
float a[n][n];
 
// Вектор правых частей
float b[n];
 
// Вектор решения
float x[n];
float xk[n],p[n];
 
 
// Выбираем начальное приближение
xk[0]=1;xk[1]=1;
 
do
{
 
iter++;
 
for (i=0;i<n;i++)
{
x[i] = xk[i];
}
 
// Заполняем матрицу A
a[0][0]=df1_dx(x); a[0][1]=df1_dy(x);
a[1][0]=df2_dx(x); a[1][1]=df2_dy(x);
 
// Заполняем вектор правых частей
b[0]=-f1(x); b[1]=-f2(x);
 
GaussSolve (a,b,p);
 
// Получаем следующее приближение,
// складывая предыдущее и вновь полученное
for (i=0;i<n;i++)
{
xk[i] = x[i] + p[i];
}
 
}
while (fabs(norm (x) - norm(xk))>=eps);
 
 
printf ("%s","Вектор решения: [ ");
for (i=0;i<n;i++)
{
printf ("%f ",xk[i]);
}
printf ("%s","]\n");
printf ("Решение найдено за %i интераций\n",iter);
 
 
return 0;
}
 
 
float f1(float x[n])
{
return x[0]*x[1]+x[1]*x[1]+2;
}
 
float f2(float x[n])
{
return x[0]*x[1]*x[1]-3;
}
 
float df1_dx(float x[n])
{
return x[1];
}
 
float df1_dy(float x[n])
{
return x[0]+2*x[1];
}
 
float df2_dx(float x[n])
{
return x[1]*x[1];
}
 
float df2_dy(float x[n])
{
return 2*x[0]*x[1];
}
 
float norm (float x[n])
{
float sum = 0;
 
for (int i=0;i<n;i++)
{
sum+=x[i]*x[i];
}
 
return sqrtf(sum);
}
 
void GaussSolve (float a[n][n],float b[n],float x[n])
 
// Решение СЛУ методом Гаусса с выбором главного элемента по столбцу
// На входе: матрица коэффициентов, вектор правых частей, вектор решения
// На выходе: заполненный вектор решения
 
{
 
const float eps=0.0001f;
 
float tmpValue;
int i,j,k,z;
float dblLeadElement;
 
 
// Прямой ход
for(i=0; i<n; i++)
{
 
dblLeadElement=a[i][i];
 
// Находим строку, в которой элемент, стоящий под ведущим - наибольший
 
float tmpMax = dblLeadElement;
int tmpMaxNumber = i;
 
for (z=i;z<n;z++)
{
if (a[z][i]>tmpMax){tmpMax = a[z][i];tmpMaxNumber=z;}
}
 
// Меняем местами i-ю строку и строку tmpMaxNumber
for (z=i;z<n;z++)
{
tmpValue = a[i][z];
a[i][z] = a[tmpMaxNumber][z];
a[tmpMaxNumber][z] = tmpValue;
 
}
tmpValue = b[i];
b[i] = b[tmpMaxNumber];
b[tmpMaxNumber] = tmpValue;
 
 
dblLeadElement = tmpMax;
 
 
 
for(j=i; j<n; j++)
{
a[i][j]/=dblLeadElement;
}
b[i]/=dblLeadElement;
 
 
 
 
for(k=i+1; k<n; k++)
{
 
float dblToDivide=a[k][i]/a[i][i];
for(z=i;z<n; z++)
{
a[k][z]-=a[i][z]*dblToDivide;
}
 
b[k]-=b[i]*dblToDivide;
}
 
 
}
 
 
// Обратный ход
x[n-1]=b[n-1];
 
for(k=n-2; k>=0; k--)
{
float sum=b[k];
 
for(j=k+1; j<n; j++)
{
sum-=a[k][j]*x[j];
}
x[k]=sum;
}
 
}

#include <cstdlib>
#include <iostream>
#include <math.h>
#include <iomanip>
 
using namespace std;
double e = 0.001; /задаем требуемую точность
 
/*ОПИСАНИЕ ФУНКЦИИ*/
double funct(double x)
{
    double f = tan(3.1 * x) + 3 * x;
    return f;
}
 
/*ПРОИЗВОДНАЯ ФУНКЦИИ*/
double diff(double x)
{
    double f = 6.1 + 3.1 * pow(tan(3.1 * x), 2);
    return f;
}
 
/*МЕТОД НЬЮТОНА*/
int task_2()
{
    int n = 0;
    double xm = 0.95;
    double xn = 0;
    double acc = fabs(xn - xm);
    
    cout << "Решение методом Ньютона:" << endl << endl;
    while (acc > e)
    {
        n++;
        xn = xm - funct(xm)/diff(xm);
        acc = fabs(xn - xm);
        cout << setw(2) << n << " " << xm << " " << xn << " " << acc << endl;
        xm = xn;
    }
    cout << endl << "Корень уравнения x = " << xn << "±" << e << endl;
}
 
int main()
{
    setlocale(LC_ALL,"Rus");
    int ch;
   
    cout.precision(4);
    cout.setf(ios::fixed);
    
    task_2();
        
    system("PAUSE");
    return EXIT_SUCCESS;
}

// Количество уравнений
const int n = 2;

float f1(float x[n])
{
return x[0]*x[1]+x[1]*x[1]+2;
}
 
float f2(float x[n])
{
return x[0]*x[1]*x[1]-3;
}
 
float df1_dx(float x[n])
{
return x[1];
}
 
float df1_dy(float x[n])
{
return x[0]+2*x[1];
}
 
float df2_dx(float x[n])
{
return x[1]*x[1];
}
 
float df2_dy(float x[n])
{
return 2*x[0]*x[1];
}

#include <iostream.h>
#include <math.h>
double Nuton(double*p,int n,double*x,double X)
{
double s=p[n-1];
for(int i=n-2;i>=0;i--) s=p[i]+(X-x[i])*s;
return s;
}
void Nuton1(double*x,double*y,int n,double*p)
{
for(int j=0;j<n-1;j++) {
p[j]=(y[j+1]-y[j])/(x[j+1]-x[j]);
// cout<<p[j]<<"\t";
}
// cout<<"\n";
int m=n-1;
int l=2;
int k;
for(int i=n-2;i>=1;i--)
{
k=0;
for(j=1;j<=i;j++)
{
p[m]=(p[m-i]-p[m-i-1])/(x[k+l]-x[k]);
// cout<<p[m]<<"\t";
m++;
k++;
}
// cout<<"\n";
l++;
}
j=0;
for(i=1;i<n;i++)
{
y[i]=p[j];
j=j+(n-i);
}
 
}
void main()
{
int const n=5,N=(n-1)*n/2;
double x[]={1,2,3,4,5};
double y[n];
for(int i=0; i<n; i++) y[i]=1+2*x[i]+3*pow(x[i],2)+4*pow(x[i],3)+5*pow(x[i],4);
double *p=new double[N];
Nuton1(x,y,n,p);
//for (i=0;i<n;i++) cout<<y[i]<<'\t';
double z=1.5;
cout<<Nuton(y,n,x,z)<<'\t'<<1+2*z+3*pow(z,2)+4*pow(z,3)+5*pow(z,4)<<'\n';
 
}