-> Presenta 2 maneras de resolver ecuaciones diferenciales. El método de Runge Kutta de 4to Orden y el método de Diferencias Finitas.
-> Tiene una precisión comprobada de 9 decimas.
-> Las ecuaciones se declaran de forma múltiple, dentro de una misma función y utilizando una función independiente. una para cada método.
-> El programa es algo complejo de entender sin embargo tiene nombres en las funciones bastante intuitivos.
-> Este condigo contiene un bloque especial que permite un escalonamiento completo o parcial en un determinante n*n.
Acá les dejo el código:
#include <iostream>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <iomanip>
#include <string.h>
using namespace std;
double e=2.718281828459045;
double pi=3.14159265358979;
///---------------------------------------------------------------------
struct puntos
{
double x,y;
};
puntos lim[100];
///---------------------------------------------------------------------
struct extras
{
double Rx[100],Qx[100],Px[100];
};
///---------------------------------------------------------------------
struct Tridente
{
double x[100],y[100],z[100];
};
///---------------------------------------------------------------------
struct matrix
{
double m[100][100];
};
///---------------------------------------------------------------------
struct Vector
{
double v[100];
};
///---------------------------------------------------------------------
void showM(int t, matrix f,Vector res)
{
printf("\n\n La matriz resultante es: \n");
for(int i=1 ; i<=t ; i++)
{
for(int j=1 ; j<=t ; j++)
{
printf(" %.10lf ",f.m[i][j]);
}
printf(" | %.10lf \n",res.v[i]);
}
}
///---------------------------------------------------------------------
void showV(int t,Vector res)
{
printf("\n\n El vector resultante es: \n");
for(int i=1 ; i<=t ; i++)
{
printf("\n %.10lf",res.v[i]);
}
}
///---------------------------------------------------------------------
void showVCompleto(int t,Vector res)
{
printf("\n\n El vector resultante es: \n");
for(int i=0 ; i<=t ; i++)
{
printf("\n %.10lf",res.v[i]);
}
}
///---------------------------------------------------------------------
void identidad(matrix &res,int sizem)
{
for(int i=1 ; i<=sizem ; i++)
{
for(int j=1 ; j<=sizem ; j++)
{
if(i==j)
res.m[i][j]=1;
else
res.m[i][j]=0;
}
}
}
///---------------------------------------------------------------------
void finale(int sizem, double h, Vector xs, extras &uwu, matrix &aux,double yo, double yf,Vector &misol)
{
Vector AB;
double p1,p2;
double x;
for(int i=1 ; i<=sizem ; i++)
{
x=xs.v[i];
/*
uwu.Rx[i]=(double) (x*x)-1;
uwu.Qx[i]=(double) -(2/uwu.Rx[i]);
uwu.Px[i]=(double) -(uwu.Qx[i]*x);
*/
/*i*/ uwu.Rx[i]=(double) x*pow(e,2*x);
/*y*/ uwu.Qx[i]=(double) -(1+(2/(x*x)));
/*y'*/ uwu.Px[i]=(double) 2/x;
for(int j=1 ; j<=sizem ; j++)
{
if(i==j)
{
aux.m[i][j]=(double) 2 + (uwu.Qx[i]*pow(h,2));
aux.m[i][j+1]=(double) -1 + ((h/2)*uwu.Px[i]);
aux.m[i][j-1]=(double) -1 - ((h/2)*uwu.Px[i]);
}
}
p1=(double) -(h*h) * uwu.Rx[i];
if(i==1)
{
p2=(double) 1 + (h*0.5*uwu.Px[i]);
misol.v[i]=(double) p1 + (p2 * yo);
}
if(i==sizem)
{
p2=(double) 1 - (h*0.5*uwu.Px[i]);
misol.v[i]=(double) p1 + (p2 * yf);
}
if(i!=1 && i!=sizem)
{
misol.v[i]=p1;
}
///printf("\n\n Los pqr son : Px= %.10f Qx= %.10f Rx= %.10f",uwu.Px[i],uwu.Qx[i],uwu.Rx[i]);
}
}
///---------------------------------------------------------------------
void escalonar(int sizem,matrix &coe,Vector &sol)
{
Vector mul;
double num;
num=coe.m[1][1];
for(int i=1 ; i<=sizem ; i++)
{
coe.m[1][i]=(double) coe.m[1][i]/num;
}
sol.v[1]=(double) sol.v[1]/num;
for(int i=1 ; i<=sizem ; i++)
{
mul.v[i]=-coe.m[i+1][i];
sol.v[i+1]=(double) sol.v[i+1] + (sol.v[1]*mul.v[i]);
for(int j=1 ; j<=sizem ; j++)
{
coe.m[i+1][j]=(double) coe.m[i+1][j] + (coe.m[i][j]*mul.v[i]);
}
}
num=coe.m[2][2];
for(int i=2 ; i<=sizem ; i++)
{
coe.m[2][i]=(double) coe.m[2][i]/num;
}
sol.v[2]=(double) sol.v[2]/num;
for(int i=2 ; i<=sizem ; i++)
{
mul.v[i]=-coe.m[i+1][i];
sol.v[i+1]=(double) sol.v[i+1] + (sol.v[1]*mul.v[i]);
for(int j=2 ; j<=sizem ; j++)
{
coe.m[i+1][j]=(double) coe.m[i+1][j] + (coe.m[i][j]*mul.v[i]);
}
}
}
///---------------------------------------------------------------------
void escalonartotal(int sizem,matrix &coe,Vector &sol)
{
Vector mul;
double num;
for(int k=1 ; k<=sizem ; k++)
{
num=coe.m[k][k];
for(int i=k ; i<=sizem ; i++)
{
coe.m[k][i]=(double) coe.m[k][i]/num;
}
sol.v[k]=(double) sol.v[k]/num;
for(int i=k ; i<=sizem ; i++)
{
mul.v[i]=-coe.m[i+1][i];
sol.v[i+1]=(double) sol.v[i+1] + (sol.v[k]*mul.v[i]);
for(int j=i ; j<=sizem ; j++)
{
coe.m[i+1][j]=(double) coe.m[i+1][j] + (coe.m[i][j]*mul.v[i]);
}
}
}
///Sacar inversa, cambiando el vector solucion para asi cambiar el resultado de manera directa.
///Sin la necesidad de reemplazar valores en forma ascendente.
for(int k=sizem ; k>=1 ; k--)
{
for(int i=k ; i>=1 ; i--)
{
mul.v[i]=-coe.m[i-1][i];
sol.v[i-1]=(double) sol.v[i-1] + (sol.v[k]*mul.v[i]);
for(int j=i ; j>=1 ; j--)
{
coe.m[i-1][j]=(double) coe.m[i-1][j] + (coe.m[i][j]*mul.v[i]);
}
}
}
}
///---------------------------------------------------------------------
void metdiffin()///METODO DIFERENCIAS FINITAS
{
system("cls");
int NI,i=0;
double a,b,c,d,h;
matrix res;
extras ext;
Vector auxi;
Vector solucion;
printf("\n Metodo de Diferencias Finitas\n");
printf("\n Ingrese numero de intervalos para resolver: "); cin>>NI;
printf("\n Coloque los puntos inicial y final \n");
//printf(" x0-> ");cin>>lim[0].x;
printf(" y0-> ");cin>>lim[0].y;
printf("\n");
//printf(" xf-> ");cin>>lim[NI].x;
printf(" yf-> ");cin>>lim[NI].y;
lim[0].x=pi/2;
lim[NI].x=pi;
a=lim[0].x;
b=lim[NI].x;
c=lim[0].y;
d=lim[NI].y;
h=(double) (b-a)/NI;
printf("\n Paso h -> %.10lf\n",h);
do
{
lim[i].x=a+(i*h);
printf("\n x[%d] = %.10lf",i,lim[i].x);
auxi.v[i]=lim[i].x;
i++;
}while(lim[i-1].x != b);
printf("\n\n y[0] = %.10lf",c);
printf("\n y[%d] = %.10lf",NI,d);
identidad(res,NI-1);
finale(NI-1,h,auxi,ext,res,c,d,solucion); ///en ext esta PQR
showM(NI-1,res,solucion);
//escalonar(NI-1,res,solucion);
escalonartotal(NI-1,res,solucion);
showM(NI-1,res,solucion);
for(int t=1 ; t<NI ; t++)
{
lim[t].y=solucion.v[t];
}
printf("\n\n");
system("PAUSE");
system("cls");
printf("\n Tabla de Resultados");
printf("\n---------------------------------------------");
for(int t=0 ; t<=NI ; t++)
{
printf("\n %d. | %.10lf | %.10lf ",t,lim[t].x,lim[t].y);
}
printf("\n\n");
}
///---------------------------------------------------------------------
void gfuncion(double x, double y , double z,double &gfun)
{
///ydd -> y'' -> z' ->g(x,y,z)
///yd -> y' -> z
///gfun= -> y''
gfun = (double) (2+(1/x))*z - (1+(1/x))*y + (x+1-(1/x))*pow(e,2*x);
//gfun = (double) (2/x*z) - ((2*y)/(x*x)) + (x*log(x));
}
///---------------------------------------------------------------------
void Trans(int sizem,double h,Tridente &Teseo)
{
double k1,k2,k3,k4,l1,l2,l3,l4;
double xx,yy,zz,gf,p1,p2,p=(double) 1/6,p3,p4;
for(int i=0 ; i<=sizem ; i++)
{
xx=Teseo.x[i];
yy=Teseo.y[i];
zz=Teseo.z[i];
k1=(double) h*zz;
gfuncion(xx , yy , zz,gf);
l1=(double) h*gf;
k2=(double) h*(zz+(l1/2));
gfuncion( xx+(h/2) , yy+(k1/2) , zz+(l1/2) ,gf );
l2=(double) h*gf;
k3=(double) h*(zz+(l2/2));
gfuncion( xx+(h/2) , yy+(k2/2) , zz+(l2/2) ,gf );
l3=(double) h*gf;
k4=(double) h*(zz+l3);
gfuncion( xx+h , yy+k3 , zz+l3 ,gf);
l4=(double) h*gf;
//printf("\n K| %.10lf | %.10lf | %.10lf | %.10lf",k1,k2,k3,k4);
//printf("\n L| %.10lf | %.10lf | %.10lf | %.10lf\n",l1,l2,l3,l4);
p1=(double) k1 + (2*k2) + (2*k3) + k4;
p2=(double) l1 + (2*l2) + (2*l3) + l4;
p3=(double) p*p1;
p4=(double) p*p2;
Teseo.y[i+1]=(double) yy + p3;
Teseo.z[i+1]=(double) zz + p4;
/*printf(" Parte sumandos\n");
printf(" K -> %.10lf",p1);
printf(" L -> %.10lf",p2);
printf(" P -> %.10lf",p);
printf(" L -> %.10lf",p2);
printf("\n\n");*/
}
}
///---------------------------------------------------------------------
void Condiciones(int sizem,Tridente &Perseo)
{
///valor inicial de y, con x inicial
double y=1.51977;
//double y=1;
///valor inicial de y, con x inicial, funcion derivada
double yd=3.30996;
//double yd=0;
Perseo.y[0]=y;
Perseo.z[0]=yd;
}
///---------------------------------------------------------------------
void showVVV(int sizem,Tridente Triton)
{
printf("\n\n n | X | Y | Z\n");
printf("-------------------------------------------------");
for(int j=0 ; j<=sizem ; j++)
{
printf("\n %d.| %.10lf | %.10lf | %.10lf",j,Triton.x[j],Triton.y[j],Triton.z[j]);
}
}
///---------------------------------------------------------------------
void metrungekutta()
{
system("cls");
int NI,i=0;
double a,b,h;
Tridente Poseidon;
printf("\n METODO DE RUNGE KUTTA\n");
printf("\n Ingrese numero de intervalos para resolver: "); cin>>NI;
printf("\n Coloque el intervalo de X \n");
printf(" xo-> "); cin>>a;
printf(" yf-> "); cin>>b;
printf("\n");
h=(double) (b-a)/NI;
printf("\n Paso h -> %.10lf\n",h);
///GENERANDO XS
do
{
Poseidon.x[i] = a + (i*h);
i++;
}while(Poseidon.x[i-1] != b);
///GENERANDO YS = 0
for(int j=0 ; j<=NI ; j++)
{
Poseidon.y[j] = 0;
}
///GENERANDO ZS = 0
for(int j=0 ; j<=NI ; j++)
{
Poseidon.z[j] = 0;
}
printf("\n\n LOS VALORES INICIALES SON:");
Condiciones(NI,Poseidon);
showVVV(NI,Poseidon);
printf("\n\n LOS VALORES FINALES SON:");
Trans(NI,h,Poseidon);
showVVV(NI,Poseidon);
}
///---------------------------------------------------------------------
void menu()
{
int des;
printf("\n Cual de los metodos desea utilizar:\n");
printf("\n 1. Metodo de Runge-Kutta de 4to Orden");
printf("\n 2. Metodo de Diferencias Finitas");
printf("\n\n -> "); cin>>des;
if(des==1) metrungekutta();
else metdiffin();
}
///---------------------------------------------------------------------
int main()
{
menu();
return 0;
}
Comments