// w.cpp : Defines the entry point for the console application.
//
#include "stdio.h"
#include "stdlib.h"
#include "math.h"
#define NUM 3

void gauss (double A[][NUM], double b[], double x[])
{
int i, j, k;
int n=NUM;

for (k=0; k {
if (A[k][k] == 0.0)
for (i=k+1; i if (A[i][k] != 0.0)
{
double tmp;
for (j=0; j {
tmp = A[k][j];
A[k][j] = A[i][j];
A[i][j] = tmp;
}
tmp = b[k];
b[k] = b[i];
b[i] = tmp;
break;
}

for (i=k+1; i {
double m = A[i][k] / A[k][k];
for (j=0; j {
A[i][j] -= m * A[k][j];
}
b[i] -= m * b[k];
}
}

for (i=n-1; i>=0; i--)
{
double sum = 0.0;
for (j=n-1; j>=i+1; j--)
sum += A[i][j]*x[j];
x[i] = (b[i] - sum)/A[i][j];
}
}
void F(double x[], double wart_F[])
{
wart_F[0]=x[0]*x[0]-x[1]*x[1]-4*x[2];
wart_F[1]=x[0]*x[0]+x[1]*x[1]-2;
wart_F[2]=-x[0]+x[1]*x[1]+x[2]*x[2];
return;
}
void Jakobian(double x[], double wart_Jakob[][NUM])
{
wart_Jakob[0][0]=2.*x[0];
wart_Jakob[0][1]=-2.*x[1];
wart_Jakob[0][2]=-4.;
wart_Jakob[1][0]=2*x[0];
wart_Jakob[1][1]=2.*x[2];
wart_Jakob[1][2]=0.;
wart_Jakob[2][0]=-1.;
wart_Jakob[2][1]=2.*x[1];
wart_Jakob[2][2]=2.*x[2];

return;
}
int main()
{
int n=3; //liczba równań
int i=1; //max liczba iteracji
double x[3], wart_F[3], wart_Jakob[3][3];
double eps, norma, sig=0.1;
double popr[3];

printf(" podaj dokladniosc rozwiazania epsilon=");
scanf("%lg",&eps);
printf(" podaj 3-elementowy wektor poczatkowy x=");
scanf("%lg", &x[0]);
scanf("%lg", &x[1]);
scanf("%lg", &x[2]);

F(x, wart_F);
norma=sqrt(wart_F[0]*wart_F[0]+wart_F[1]*wart_F[1]+wart_F[2]*wart_F[2]);
while(norma>eps && i<1000 )
{
Jakobian(x, wart_Jakob);
wart_F[0]=-1*wart_F[0];
wart_F[1]=-1*wart_F[1];
wart_F[2]=-1*wart_F[2];


gauss(wart_Jakob, wart_F, popr);

x[0]=x[0]+sig*popr[0];
x[1]=x[1]+sig*popr[1];
x[2]=x[2]+sig*popr[2];

F(x, wart_F);
norma=sqrt(wart_F[0]*wart_F[0]+wart_F[1]*wart_F[1]+wart_F[2]*wart_F[2]);
i++;
}
printf("Ilosc iteracji: %d\n", i);
printf("Rozwiazanie ukladu:\n %lg\n %lg\n %lg\n", x[0], x[1], x[2]);
system("PAUSE");
return 0;
}