C++ GAUSS J. Notas de Clase.

NOTAS DE CLASE.

Programa realizado en C++ durante la clase, Análisis Numérico por el profesor Hector Hernadez.

#include<math.h>

#include<stdio.h>

#include<stdlib.h>

#define N 100

#define Max_np 1000 

void ImpMat(double M[N][N+1],int n, char Let[10]){

int i,j;

printf("\n%s\n",Let);

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

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

printf("%lf ",M[i][j]);

}

printf("\n");

}

}

int LugarMA(double M[N][N+1],int n,int piv){

int lugar=piv,i;

double epsilon=0.000001;

double MayVA=fabs(M[piv][piv]);

for(i=piv+1;i<n;i++){

if(fabs(M[i][piv])>MayVA){

MayVA=fabs(M[i][piv]);

lugar=i;

}

}

if(MayVA<epsilon) lugar=-1;

return lugar;

}

void InterRen(double M[N][N+1],int n,int r1, int r2){

int j;

double temp;

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

  temp=M[r1][j];

  M[r1][j]=M[r2][j];

  M[r2][j]=temp;

}

  

}

int Pivotaje(double M[N][N+1],int n,int piv){

double temp;

int i,j,Renglon;

Renglon=LugarMA(M,n,piv);

//printf("\n EN contré el renglon %d\n",Renglon);

if(Renglon>=0){

    if(Renglon > piv) InterRen(M,n,Renglon,piv);

            temp=1/M[piv][piv];

//Hace el uno

for(j=0;j<=n;j++){                 // j++      j=j+1

M[piv][j]=M[piv][j]*temp;

}  

//Hace los ceros de arriba del uno  

for(i=0;i<piv;i++){

    temp=M[i][piv];

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

    M[i][j]=M[i][j]-M[piv][j]*temp;

}    

}   

    //hace los ceros de abajo del uno

for(i=piv+1;i<n;i++){

    temp=M[i][piv];

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

    M[i][j]=M[i][j]-M[piv][j]*temp;

}    

} 

return 0;

}else{

return -1;

}

 

}

void GaussJ(double M[N][N+1],int n,int *SU){

int piv;

for(piv=0;piv<n;piv++){ //Cuando quiero intentar continuar con pivotajes

if(Pivotaje(M,n,piv)<0) *SU=0;

//ImpMat(A,orden,"A=");

}

}

int main(){

    double J[N][N+1];

    int orden=3;

double x[Max_np],y[Max_np],a,b,c,xic,xi2c,lnxi,da,db,dc;

int i,np=5;

int SolU=1,k;

  x[0]=.2 ; y[0]=8.03720209;

x[1]=.7 ; y[1]=8.28524737;

x[2]=1 ; y[2]=8.88847557;

x[3]=2 ; y[3]=17.2462463;

x[4]=2.7 ; y[4]=33.8227762;

a=5.91;

b=3.87;

c=1.94;

J[0][0]=np; J[0][1]=0; J[0][2]=0; J[0][3]=0;

J[1][0]=0;  J[1][1]=0; J[1][2]=0; J[1][3]=0;

J[2][0]=0;  J[2][1]=0; J[2][2]=0; J[2][3]=0;

for(i=0;i<np;i++){

xic=pow(x[i],c);

lnxi=log(x[i]);

xi2c=pow(x[i],2*c);

J[0][1]=J[0][1]+xic;

J[0][2]=J[0][2]+b*xic*lnxi;

J[0][3]=J[0][3]+a+b*xic-y[i];

J[1][0]=J[1][0]+xic;

J[1][1]=J[1][1]+xi2c;

J[1][2]=J[1][2]+(a+2*b*xic-y[i])*xic*lnxi;

J[1][3]=J[1][3]+a*xic+b*xi2c-xic*y[i];

J[2][0]=J[2][0]+xic*lnxi;

J[2][1]=J[2][1]+xi2c*lnxi;

J[2][2]=J[2][2]+b*xi2c*2*lnxi*lnxi;

J[2][3]=J[2][3]+a*xic*lnxi+b*xi2c*lnxi-xic*y[i]*lnxi;

}

ImpMat(J,orden,"SEL=");

system("pause");

//Newton

for(k=0;k<3000;k++){

GaussJ(J,orden,&SolU);

da=J[0][3];

db=J[1][3];

dc=J[2][3];

printf("\n deltas %lf, %lf ,%lf",da,db,dc);

a=a-da;

b=b-db;

c=c-dc;

J[0][0]=np; J[0][1]=0; J[0][2]=0; J[0][3]=0;

J[1][0]=0;  J[1][1]=0; J[1][2]=0; J[1][3]=0;

J[2][0]=0;  J[2][1]=0; J[2][2]=0; J[2][3]=0;

for(i=0;i<np;i++){

xic=pow(x[i],c);

lnxi=log(x[i]);

xi2c=pow(x[i],2*c);

J[0][1]=J[0][1]+xic;

J[0][2]=J[0][2]+b*xic*lnxi;

J[0][3]=J[0][3]+a+b*xic-y[i];

J[1][0]=J[1][0]+xic;

J[1][1]=J[1][1]+xi2c;

J[1][2]=J[1][2]+(a+2*b*xic-y[i])*xic*lnxi;

J[1][3]=J[1][3]+a*xic+b*xi2c-xic*y[i];

J[2][0]=J[2][0]+xic*lnxi;

J[2][1]=J[2][1]+xi2c*lnxi;

J[2][2]=J[2][2]+b*xi2c*2*lnxi*lnxi;

J[2][3]=J[2][3]+a*xic*lnxi+b*xi2c*lnxi-xic*y[i]*lnxi;

}

ImpMat(J,orden,"SEL=");

printf("\n |F|=%lf",sqrt(J[0][3]*J[0][3]+J[1][3]*J[1][3]+J[2][3]*J[2][3]));

}

printf("\n a=%lf b=%lf c=%lf",a,b,c);