-> Mantiene una precisión comprobada de 9 decimas.
-> Presenta 1 métodos para resolver el sistema planteado, el método de Newton-Rhapson.
-> Las ecuaciones y derivadas se trabajan en funciones independientes.
-> Las ecuaciones y derivadas deben ser escritas antes de correr el programa.
-> Hay una variante dentro del código mejorable, el trabajo con matrices como escalonamiento. En el código de Resolver Ecuaciones Diferenciales esta la variante superior que permite trabajar un determinante n*n
-> El programa también muestra una solución detallada con las matrices cambiadas.
-> El programa presenta al final una tabla de valores y su error durante sus iteraciones.
-> Hará cuantas iteraciones sean necesarias para obtener un error menor a 0.000000000
Acá les dejo el código:
#include <iostream>
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
using namespace std;
double pi=3.14159265358979;
double e=2.718281828459045;
///---------------------------------------------------------------------
struct Matriz
{
double m[2][2];
};
struct Vect
{
double v[10];
};
struct Cadena
{
double x[10],y[10];
};
///---------------------------------------------------------------------
void Identidad(int lm , Matriz &a)
{
for(int i=0 ; i<lm ; i++)
{
for(int j=0 ; j<lm ; j++)
{
if(i==j)
a.m[i][j]=1;
else
a.m[i][j]=0;
}
}
}
///---------------------------------------------------------------------
double miabs(double &val)
{
if(val<0)
val=val*(-1);
return val;
}
///---------------------------------------------------------------------
void derivadas(double &f1x, double &f2x, double &f1y, double &f2y, double x, double y) //x0 = 1; y0 = 1 ; Tol=0.0000
{
f1x=(double) -3 * pow((x-1),2); f1y=(double) 2*y;
f2x=(double) -pow(e,-x) - 3; f2y=(double) ((2*y) / (1+(y*y))) + 2;
}
///---------------------------------------------------------------------
void ecuaciones(double &f1, double &f2, double x, double y)
{
//f1=(double) pow((x+2),2) + (4*y*y) - 4;
//f2=(double) (4*x*x) - pow((y-1),2) - 4;
f1=(double) pow(y,2) - pow((x-1),3);
f2=(double) log(1+pow(y,2)) + pow(e,-x) - (3*x) + (2*y);
}
///---------------------------------------------------------------------
void vm(Matriz &res)
{
printf("\n\n\n\n");
for(int i=0 ; i<2 ; i++)
{
for(int j=0 ; j<2 ; j++)
{
printf(" %.7f ",res.m[i][j]);
}
cout<<endl;
}
}
///---------------------------------------------------------------------
void inversa(Matriz &c,Matriz &ide)
{
double vali=c.m[0][0],vali2=c.m[1][0];
for(int i=0 ; i<2 ; i++)
{
c.m[0][i]=(double) c.m[0][i]/vali;
ide.m[0][i]=(double) ide.m[0][i]/vali;
c.m[1][i]=(double) c.m[1][i]/vali2;
ide.m[1][i]=(double) ide.m[1][i]/vali2;
}
for(int i=0 ; i<2 ; i++)
{
c.m[1][i]=(double) c.m[1][i]+(-1*c.m[0][i]);
ide.m[1][i]=(double) ide.m[1][i]+(-1*ide.m[0][i]);
}
double vali3=c.m[1][1];
for(int i=0 ; i<2 ; i++)
{
c.m[1][i]=(double) c.m[1][i]/vali3;
ide.m[1][i]=(double) ide.m[1][i]/vali3;
}
double vali4=c.m[0][1];
for(int i=0 ; i<2 ; i++)
{
c.m[0][i]=(double) c.m[0][i]-(c.m[1][i]*vali4);
ide.m[0][i]=(double) ide.m[0][i]-(ide.m[1][i]*vali4);
}
}
///---------------------------------------------------------------------
void inversa2(Matriz &num,Matriz &ide)
{
double vali = num.m[0][0];
for(int i=0 ; i<2 ; i++)
{
num.m[0][i]=(double) num.m[0][i] / vali;
ide.m[0][i]=(double) ide.m[0][i] / vali;
}
double vali2 = -num.m[1][0];
for(int j=0 ; j<2 ; j++)
{
num.m[1][j]=(double) num.m[1][j] + (vali2*num.m[0][j]);
ide.m[1][j]=(double) ide.m[1][j] + (vali2*ide.m[0][j]);
}
double vali3 = num.m[1][1];
for(int k=0 ; k<2 ; k++)
{
num.m[1][k]=(double) num.m[1][k] / vali3;
ide.m[1][k]=(double) ide.m[1][k] / vali3;
}
double vali4 = -num.m[0][1];
for(int l=0 ; l<2 ; l++)
{
num.m[0][l]=(double) num.m[0][l] + (vali4*num.m[1][l]);
ide.m[0][l]=(double) ide.m[0][l] + (vali4*ide.m[1][l]);
}
}
///---------------------------------------------------------------------
void matxvect(Matriz &res, Vect &fs)
{
double aux0=fs.v[0],aux1=fs.v[1];
fs.v[0]=(double) (res.m[0][0]*aux0)+(res.m[0][1]*aux1);
fs.v[1]=(double) (res.m[1][0]*aux0)+(res.m[1][1]*aux1);
}
///---------------------------------------------------------------------
void vectxvect(Vect &inis, Vect &fs)
{
fs.v[0]=(double) inis.v[0]-fs.v[0];
fs.v[1]=(double) inis.v[1]-fs.v[1];
}
///---------------------------------------------------------------------
void MetNewRap()
{
int des,j=0,i=0;
Matriz coef;
Matriz ident;
Matriz respaldo;
Vect iniciales;
Vect efes;
Vect E;
Cadena aux;
Cadena xsys;
Identidad(2,ident);
Identidad(2,respaldo);
double Tol,x0,y0,f1,f2,f1x,f2x,f1y,f2y,a,b,c,d,a1,a2;
printf("\n Metodo de Newton Rhapson \n");
printf("\n Introduzca la Toleracia -> "); scanf("%lf",&Tol);
printf("\n Introduza valor x0 -> "); scanf("%lf",&x0);
printf("\n Introduza valor y0 -> "); scanf("%lf",&y0);
derivadas(f1x,f2x,f1y,f2y,x0,y0);
ecuaciones(f1,f2,x0,y0);
a=f1x; b=f1y;
c=f2x; d=f2y;
printf("\n La Matriz Jacobiana del sistema es: ");
printf("\n |%.6f %.6f| \n |%.6f %.6f|",miabs(a),miabs(b),miabs(c),miabs(d));
printf("\n\n Se debe de Reorganizar?\n 1. SI (Si A[0][0] = 0)\n 2. NO");
printf("\n -> "); cin>>des;
if(des==1)
{
printf("\n La Matriz Jacobiana reorganizada del sistema es: ");
printf("\n |%.6f %.6f| \n |%.6f %.6f|\n\n",miabs(c),miabs(d),miabs(a),miabs(b));
printf("\n Luego: \n");
printf("\n |x| = |x0| - |f2x f2y|^(-1) |f2|");
printf("\n |y| = |y0| |f1x f1y| |f1|\n\n");
system("PAUSE");
system("cls");
do
{
coef.m[0][0]=f2x; coef.m[0][1]=f2y;
coef.m[1][0]=f1x; coef.m[1][1]=f1y;
printf("\n %d. ITERACION\n",i+1);
printf("\n |x| = |%.6f| - |%.10f %.10f|^(-1) |%.10f|",x0,f2x,f2y,f2);
printf("\n |y| = |%.6f| |%.10f %.10f| |%.10f|\n",y0,f1x,f1y,f1);
inversa2(coef,ident);
printf("\n |x| = |%.6f| - |%.10f %.10f| |%.10f|",x0,ident.m[0][0],ident.m[0][1],f2);
printf("\n |y| = |%.6f| |%.10f %.10f| |%.10f|",y0,ident.m[1][0],ident.m[1][1],f1);
iniciales.v[0]=x0; efes.v[0]=f2;
iniciales.v[1]=y0; efes.v[1]=f1;
matxvect(ident,efes);
ident=respaldo;
vectxvect(iniciales,efes);
printf("\n\n RESULTADOS");
printf("\n |x| = %.10lf",efes.v[0]);
printf("\n |y| = %.10lf",efes.v[1]);
xsys.x[i]=efes.v[0];
xsys.y[i]=efes.v[1];
ecuaciones(a1,a2,xsys.x[i],xsys.y[i]);
E.v[i]=(float) miabs(a1) + miabs(a2);
printf("\n |E| = %.10lf\n",E.v[i]);
x0=xsys.x[i];
y0=xsys.y[i];
derivadas(f1x,f2x,f1y,f2y,x0,y0);
ecuaciones(f1,f2,x0,y0);
i++;
}while(Tol<E.v[i-1]);
system("PAUSE");
system("cls");
printf("\n TABLA DE RESULTADOS\n");
printf(" -------------------------------------------------");
printf("\n n x y E\n");
printf(" -------------------------------------------------");
while(E.v[j]> Tol)
{
printf("\n %d. | %.10lf | %.10lf | %.10lf",j+1,xsys.x[j],xsys.y[j],E.v[j]);
j++;
}
printf("\n %d. | %.10lf | %.10lf | %.10lf",j+1,xsys.x[j],xsys.y[j],E.v[j]);
cout<<"\n"<<endl;
system("PAUSE");
}
if(des==2)
{
printf("\n Luego: \n");
printf("\n |x| = |x0| - |f1x f1y|^(-1) |f1|");
printf("\n |y| = |y0| |f2x f2y| |f2|\n\n");
system("PAUSE");
system("cls");
do
{
coef.m[0][0]=f1x; coef.m[0][1]=f1y;
coef.m[1][0]=f2x; coef.m[1][1]=f2y;
printf("\n %d. ITERACION\n",i+1);
printf("\n |x| = |%.6f| - |%.10f %.10f|^(-1) |%.10f|",x0,f1x,f1y,f1);
printf("\n |y| = |%.6f| |%.10f %.10f| |%.10f|\n",y0,f2x,f2y,f2);
inversa(coef,ident);
printf("\n |x| = |%.6f| - |%.10f %.10f| |%.10f|",x0,ident.m[0][0],ident.m[0][1],f1);
printf("\n |y| = |%.6f| |%.10f %.10f| |%.10f|",y0,ident.m[1][0],ident.m[1][1],f2);
iniciales.v[0]=x0; efes.v[0]=f1;
iniciales.v[1]=y0; efes.v[1]=f2;
matxvect(ident,efes);
ident=respaldo;
vectxvect(iniciales,efes);
printf("\n\n RESULTADOS");
printf("\n |x| = %.10lf",efes.v[0]);
printf("\n |y| = %.10lf",efes.v[1]);
xsys.x[i]=efes.v[0];
xsys.y[i]=efes.v[1];
ecuaciones(a1,a2,xsys.x[i],xsys.y[i]);
E.v[i]=(float) miabs(a1) + miabs(a2);
printf("\n |E| = %.10lf\n",E.v[i]);
x0=xsys.x[i];
y0=xsys.y[i];
derivadas(f1x,f2x,f1y,f2y,x0,y0);
ecuaciones(f1,f2,x0,y0);
i++;
}while(Tol<E.v[i-1]);
system("PAUSE");
system("cls");
printf("\n TABLA DE RESULTADOS\n");
printf(" -------------------------------------------------");
printf("\n n x y E\n");
printf(" -------------------------------------------------");
while(E.v[j]> Tol)
{
printf("\n %d. | %.10lf | %.10lf | %.10lf",j+1,xsys.x[j],xsys.y[j],E.v[j]);
j++;
}
printf("\n %d. | %.10lf | %.10lf | %.10lf",j+1,xsys.x[j],xsys.y[j],E.v[j]);
cout<<"\n"<<endl;
system("PAUSE");
}
}
///---------------------------------------------------------------------
int main()
{
MetNewRap();
return 0;
}
コメント