Universidad Nacional de Salta
Facultad de Ingeniería
Ingeniería Química
Análisis Numérico
Trabajo Practico N°8
ECUACIONES DIFERENCIALES ORDINARIAS. METODOS DE RUNGE-KUTTA
ARAMAYO, Juan Esteban.
L.U: 314.197.
PROBLEMA 8.1
Desarrolle un programa en C++ de un método de Runge-Kutta genérico de acuerdo con el seudo código listado a
continuación. Lea los datos desde un archivo y efectúe una salida de resultados en forma tabular, también en un archivo
que contenga los valores de la iteración n, el punto x y el valor de la solución aproximada Y. Este último valor se debe
mostrar con 8 dígitos decimales.
Lista de símbolos utilizados en el algoritmo
• k : orden del método • N : número de puntos en donde se quiere determinar la solución
• h : longitud del paso • (x0, Y0 ) : valores iniciales
• ai, bi, ci,j: coeficientes del método • S : vector auxiliar
• x, Y: vectores • n : variable auxiliar
• R : vector auxiliar
/*TRABAJO PRACTICO 8*/
/*Ecuaciones Diferenciales Ordinarias de Primer Orden. Metodos de Runge-Kutta*/
#include"..\anqlib++\anqlib64++.h"
#include"..\errores.h"
//DECLARACION DE PROTOTIPOS DE FUNCIONES
int Obtener_Datos(void);
int Aplicar_Metodo(void);
void MuestraResultados_en_Ventana(int error);
int Salida_en_Archivo(int error);
double f( double x, double y);
1
Universidad Nacional de Salta
Facultad de Ingeniería
Ingeniería Química
Análisis Numérico
// Funcion f(x,y)
double f(double x, double y)
{
return x*x*y -y; //EJERCICIO 8.1 - 8.2
return sqrt(1 + sin(x))*exp(1-y*y); //EJERCICIO 8.3
return 0.659*y*(1-y/99.0); //EJERCICIO 8.4
}
//DECLARACION DE VARIABLES GLOBALES
double a[10], b[10], c[10][10], x[100],Y[100],h ;
int k, N;
//FUNCION PRINCIPAL DEL PROGRAMA
int main()
{
int error;
//Obtener los datos del problema
error=Obtener_Datos();
//Fin datos
//Efecturar calculos
if(error==ERROR_OK) error=Aplicar_Metodo();
//Fin efecturar calculos
//Salida de resultados en archivo
if(error==ERROR_OK) error=Salida_en_Archivo(error);
//Fin salida en archivo
//Mensajes en ventana de ejecucion
MuestraResultados_en_Ventana(error);
//Fin Mensajes en ventana
2
Universidad Nacional de Salta
Facultad de Ingeniería
Ingeniería Química
Análisis Numérico
system("PAUSE");
return error;
}
//FUNCION OBTENER DATOS
int Obtener_Datos()
{
//Declracion de variables locales
int i, j, error;
char buffer_lectura[5000]; //Buffer almaena datos del archivo
//Apertura y lectura de datos
error=LeeDatos_de_Archivo("[Link]", buffer_lectura);
if(error!=0) return error;
//Obtener orden del Metodo
error=leeEntero(buffer_lectura, &k);
if(error!=0) return error;
//Obtener coeficientes a[i]
for(i=1 ; i<=k ; i++)
{
error=leeDouble(buffer_lectura ,&a[i]);
if(error!=0) return error;
}
//Obtener coeficientes b[i]
for(i=1 ; i<=k ; i++)
{
3
Universidad Nacional de Salta
Facultad de Ingeniería
Ingeniería Química
Análisis Numérico
error=leeDouble(buffer_lectura, &b[i]);
}
//Obtener coeficientes c[i][j]
for( i=2 ; i<=k ; i++)
{
for(j=1 ; j<=i-1 ; j++)
{
error=leeDouble(buffer_lectura , &c[i][j]);
if(error!=0) return error;
}
}
//Obtener cantidad N de puntos
error=leeEntero(buffer_lectura, &N);
if(error!=0) return error;
//Obtener valores iniciales
error=leeDouble(buffer_lectura, &x[0]) ;
if(error!=0) return error;
error=leeDouble(buffer_lectura, &Y[0]) ;
if(error!=0) return error;
//Obtener valor de paso h
error=leeDouble(buffer_lectura, &h) ;
if(error!=0) return error;
return ERROR_OK;
}
4
Universidad Nacional de Salta
Facultad de Ingeniería
Ingeniería Química
Análisis Numérico
int Aplicar_Metodo()
{
//Declaracion de variables locales
int i, j,n;
double S[3], R[5];
for(n=0 ; n<=N-1 ; n++)
{
S[1]=0.0;
for(i=1 ; i<=k ; i++)
{
S[2]=0.0;
for(j=1 ; j<=i-1 ; j++)
{
S[2]=S[2] +c[i][j]*R[j];
}
R[i]=f(x[n]+b[i]*h, Y[n] + h*S[2]);
S[1]=S[1] + a[i]*R[i];
}
x[n+1]= x[n] + h;
Y[n+1]= Y[n] + h*S[1];
}
return ERROR_OK;
}
//FUNCION MUESTRA RESULTADOS EN PANTALLA
void MuestraResultados_en_Ventana(int error)
{
5
Universidad Nacional de Salta
Facultad de Ingeniería
Ingeniería Química
Análisis Numérico
int error2;
if(error!=ERROR_OK)
{
//Verificar si se produjeron errores
//Errores de la librería anqlib++
error2=ImprimeTextoErrorV(error);
if(error2<0)
{
//Errores de cada método numérico en particular
if(error==ERROR_PIVOTE_PCERO){ cout<<"Error. Pivote proximo a cero.\r\n";
return; }
if(error==ERROR_NO_SOLUCION ) { cout<<"Error. Sistema sin solucion.\r\n" ;
return; }
if(error==ERROR_DERIVADA_PCERO){cout<<"Error. Derivada proxima a cero.\r\n" ;
return;}
cout<<"Error desconocido.\r\n"; return;
}
}
else
{
cout<<"Ejecucion sin errores\r\n"; //mensaje de ejecución satisfactoria
}
}
//FUNCION SALIDA EN ARCHIVO
int Salida_en_Archivo(int error)
{
int n;
char buffer_salida [50000],buff_numero[100];
6
Universidad Nacional de Salta
Facultad de Ingeniería
Ingeniería Química
Análisis Numérico
lstrcpy(buffer_salida, "n x[n] Y[n]\r\n");
for(n=0 ; n<=N ; n++)
{
ConvertirEntero_a_Texto(n, buff_numero); lstrcat(buffer_salida,
buff_numero);
lstrcat(buffer_salida, "\t");
ConvertirDouble_a_Texto(x[n],8, buff_numero); lstrcat(buffer_salida,
buff_numero);
lstrcat(buffer_salida, "\t");
ConvertirDouble_a_Texto(Y[n], 8, buff_numero); lstrcat(buffer_salida,
buff_numero);
lstrcat(buffer_salida, "\r\n");
}
//Copiar el contenido de buffer_salida en arhivo de texto de salida
error=EscribeResultados_en_Archivo("[Link]", buffer_salida);
return error;
}
Archivo de datos
7
Universidad Nacional de Salta
Facultad de Ingeniería
Ingeniería Química
Análisis Numérico
Archivo de salida
FORMA ANALITICA
𝒚′ (𝒙) = 𝒙𝟐 𝒚(𝒙) − 𝒚(𝒙) 𝒚(𝟎) = 𝟏. 𝟓
𝒅𝒚 𝒅𝒚 𝒅𝒚
= 𝒙𝟐 𝒚 − 𝒚 → = 𝒚(𝒙𝟐 − 𝟏) → = (𝒙𝟐 − 𝟏)𝒅𝒙
𝒅𝒙 𝒅𝒙 𝒚
𝒅𝒚 𝒅𝒚
∫ = ∫(𝒙𝟐 − 𝟏)𝒅𝒙 → ∫ = ∫ 𝒙𝟐 𝒅𝒙 − 𝟏 ∫ 𝒅𝒙
𝒚 𝒚
𝒙𝟑 𝒙𝟑
𝒙𝟑 ( −𝒙+𝑪) ( −𝒙+𝑪)
𝑳𝒏|𝒚| = 𝟑
− 𝒙 + 𝑪 → 𝒆𝑳𝒏|𝒚| = 𝒆 𝟑 →𝒚=𝒆 𝟑
𝟎𝟑
( −𝟎+𝑪)
𝟑
𝑺𝒐𝒍𝒖𝒄𝒊𝒐𝒏 𝑷𝒂𝒓𝒕𝒊𝒄𝒖𝒍𝒂𝒓 𝒄𝒐𝒏 𝒚(𝟎) = 𝟏. 𝟓 → 𝟏. 𝟓 = 𝒆 → 𝟏. 𝟓 = 𝒆𝑪 → 𝑳𝒏(𝟏. 𝟓) = 𝑳𝒏(𝒆𝑪 ) → 𝑪 = 𝑳𝒏(𝟏. 𝟓)
𝒙𝟑
( 𝟑 −𝒙+𝒍𝒏(𝟏.𝟓))
𝒚=𝒆
8
Universidad Nacional de Salta
Facultad de Ingeniería
Ingeniería Química
Análisis Numérico
Comparación de la solución analítica y la solución numérica en 11 puntos.
9
Universidad Nacional de Salta
Facultad de Ingeniería
Ingeniería Química
Análisis Numérico
PROBLEMA 8.2
Obtenga la solución aproximada por medio del método de Runge-Kutta de 4to orden del problema anterior para el
conjunto de coeficientes:
𝒃𝟏 = 𝟎 , 𝒃𝟐 = 𝟏/𝟑 , 𝒃𝟑 = 𝟐/𝟑 , 𝒃𝟒 = 𝟏
𝒂𝟏 = 𝒂𝟒 = 𝟏/𝟖 , 𝒂𝟐 = 𝒂𝟑 = 𝟑/𝟖
𝟏 𝟏
𝒄𝟐,𝟏 = , 𝒄𝟑,𝟏 = − , 𝒄𝟑,𝟐 = 𝟏 , 𝒄𝟒,𝟏 = 𝟏, 𝒄𝟒,𝟐 = −𝟏, 𝒄𝟒,𝟑 = 𝟏
𝟑 𝟑
Archivo de datos
Archivo de salida
10
Universidad Nacional de Salta
Facultad de Ingeniería
Ingeniería Química
Análisis Numérico
Comparacion de la solucion analitica y la numerica en 11 puntos.
11
Universidad Nacional de Salta
Facultad de Ingeniería
Ingeniería Química
Análisis Numérico
PROBLEMA 8.3
Dado el siguiente problema de valor inicial:
𝑑𝑦 2
= √1 + 𝑠𝑒𝑛(𝑥) 𝑒 (1−𝑦 )
𝑑𝑥
Para los siguientes datos:𝑥0 = 0.5, 𝑦(𝑥0 ) = −1.34, ℎ = 0.05 𝑦 𝑁 = 20, obtenga la solución aproximada por medio de los
métodos de Euler mejorado y Euler modificado.
METODO DE EULER MEJORADO
Archivo de datos
Archivo de salida
12
Universidad Nacional de Salta
Facultad de Ingeniería
Ingeniería Química
Análisis Numérico
METODO DE EULER MODIFICADO
Archivo de datos
Archivo de salida
13
Universidad Nacional de Salta
Facultad de Ingeniería
Ingeniería Química
Análisis Numérico
𝒅𝒚 𝟐
= √𝟏 + 𝒔𝒆𝒏(𝒙) 𝒆(𝟏−𝒚 )
𝒅𝒙
Separando variables se llega a:
𝒅𝒚
∫ 𝟐 = ∫ √𝟏 + 𝒔𝒆𝒏(𝒙)𝒅𝒙
𝒆𝟏−𝒚
La integral:
𝒅𝒚
∫ 𝟐
𝒆𝟏−𝒚
Es no elemental debido a esto no se puede obtener una expresión de una primitiva.
14
Universidad Nacional de Salta
Facultad de Ingeniería
Ingeniería Química
Análisis Numérico
PROBLEMA 8.4
La dinámica de reproducción de un determinado tipo de célula está dada por la siguiente ecuación diferencial:
𝒅𝒚 𝒚
= 𝟎. 𝟔𝟓𝟗𝒚(𝟏 − )
𝒅𝒙 𝟗𝟗. 𝟎
Si inicialmente se dispone de una concentración de 10 células, ¿cuántas horas serán necesarias para disponer de una
concentración de 90?. Resuelva el problema por medio del método de Runge-Kutta de 4to orden con los coeficientes
dados en la expresión (8.1) y utiliceℎ = 0.1[ℎ𝑠].
Archivo de datos
Partimos de la condición inicial de problema: para un tiempo T=0, disponemos de una concentración de células de 10
unidades.
Para obtener una concentración de 90 células, según el modelaje de la ecuación diferencial debe transcurrir un tiempo
de 6.9 horas, transcurrido ese tiempo se obtendrá una concentración de 90 células.
15
Universidad Nacional de Salta
Facultad de Ingeniería
Ingeniería Química
Análisis Numérico
Archivo de salida
16
Universidad Nacional de Salta
Facultad de Ingeniería
Ingeniería Química
Análisis Numérico
PROBLEMA 8.5
Resuelva el problema de valor inicial 8.4 mediante el método de Heun en forma manual. Utilice como datos: ℎ = 0.02,
𝑁 = 3.
El P.V.I a resolver es:
𝑑𝑦 𝑦
= 0.659𝑦(1 − )
{ 𝑑𝑡 99.0
𝑦(0) = 10
Paso ℎ = 0.02
Cantidad de puntos 𝑁 =
METODO DE HEUN
Para mayor comodidad se trabajará con la variable x=t:
𝑑𝑦 𝑦
= 0.659𝑦(1 − )
{𝑑𝑥 99.0
𝑦(0) = 10
𝟏 𝟑 𝟐 𝟐
𝒀𝒏+𝟏 = 𝒀𝒏 + 𝒉 [ 𝒇(𝒙𝒏 , 𝒀𝒏 ) + 𝒇 (𝒙𝒏 + 𝒉, 𝒀𝒏 + 𝒉. 𝒇(𝒙𝒏 , 𝒀𝒏 )) ]
𝟒 𝟒 𝟑 𝟑
𝒏 = 𝟎, 𝟏, 𝟐, … , 𝑵 − 𝟏,
Para 𝑛 = 0
𝟏 𝟑 𝟐 𝟐
𝒀𝟏 = 𝒀𝟎 + 𝒉 [ 𝒇(𝒙𝟎 , 𝒀𝟎 ) + 𝒇 ( 𝒙𝟎 + 𝒉, 𝒀𝟎 + 𝒉. 𝒇(𝒙𝟎 , 𝒀𝟎 ))]
𝟒 𝟒 𝟑 𝟑
𝟏 𝟏𝟎 𝟑 𝟐 𝟐
𝒀𝟏 = 𝟏𝟎 + 𝟎. 𝟎𝟐 [ (𝟎. 𝟔𝟓𝟗𝒙𝟏𝟎𝒙 (𝟏 − ) + 𝒇(𝟎 + 𝒙(𝟎. 𝟎𝟐), 𝟏𝟎 + 𝒙(𝟎. 𝟎𝟐)𝒙𝒇(𝟎, 𝟏𝟎)]
𝟒 𝟗𝟗 𝟒 𝟑 𝟑
𝟏 𝟏𝟎 𝟑
𝒀𝟏 = 𝟏𝟎 + 𝟎. 𝟎𝟐 [ (𝟎. 𝟔𝟓𝟗𝒙𝟏𝟎𝒙 (𝟏 − ) + 𝒇(𝟎. 𝟎𝟏𝟑𝟑𝟑𝟑, 𝟏𝟎. 𝟎𝟕𝟖𝟗𝟗𝟏 )]
𝟒 𝟗𝟗 𝟒
𝒀𝟎
𝒇(𝒙𝟎 , 𝒀𝟎 ) = 𝟎. 𝟔𝟓𝟗𝒀𝟎 (𝟏 − )
𝟗𝟗. 𝟎
17
Universidad Nacional de Salta
Facultad de Ingeniería
Ingeniería Química
Análisis Numérico
𝟏 𝟑 𝟏𝟎. 𝟎𝟕𝟖𝟗𝟗𝟏
𝒀𝟏 = 𝟏𝟎 + 𝟎. 𝟎𝟐 [ (𝟓. 𝟗𝟐𝟒𝟑𝟒𝟑) + 𝒙(𝟎. 𝟔𝟓𝟗𝒙𝟏𝟎. 𝟎𝟕𝟖𝟗𝟗𝟏𝒙(𝟏 − )]
𝟒 𝟒 𝟗𝟗. 𝟎
𝟏 𝟑
𝒀𝟏 = 𝟏𝟎 + 𝟎. 𝟎𝟐 [ (𝟓. 𝟗𝟐𝟒𝟑𝟒𝟑) + 𝒙(𝟓. 𝟗𝟔𝟓𝟖𝟒𝟎)]
𝟒 𝟒
𝒀𝟏 = 𝟏𝟎. 𝟏𝟏𝟗𝟏𝟎𝟗
Para 𝑛 = 1
𝟏 𝟑 𝟐 𝟐
𝒀𝟐 = 𝒀𝟏 + 𝒉 [ 𝒇(𝒙𝟏 , 𝒀𝟏 ) + 𝒇 ( 𝒙𝟏 + 𝒉, 𝒀𝟎 + 𝒉. 𝒇(𝒙𝟏 , 𝒀𝟏 ))]
𝟒 𝟒 𝟑 𝟑
𝟏 𝟏𝟎. 𝟏𝟏𝟗𝟏𝟎𝟗 𝟑 𝟐 𝟐
𝒀𝟏 = 𝟏𝟎. 𝟏𝟏𝟗𝟏𝟎𝟗 + 𝟎. 𝟎𝟐 [ (𝟎. 𝟔𝟓𝟗𝒙𝟏𝟎. 𝟏𝟏𝟗𝟏𝟎𝟗𝒙 (𝟏 − ) + 𝒇(𝟎 + 𝒙(𝟎. 𝟎𝟐), 𝟏𝟎. 𝟏𝟏𝟗𝟏𝟎𝟗 + 𝒙(𝟎. 𝟎𝟐)𝒙𝒇(𝟎, 𝟏𝟎. 𝟏𝟏𝟗𝟏𝟎𝟗)]
𝟒 𝟗𝟗 𝟒 𝟑 𝟑
𝟏 𝟏𝟎. 𝟏𝟏𝟗𝟏𝟎𝟗 𝟑
𝒀𝟏 = 𝟏𝟎. 𝟏𝟏𝟗𝟏𝟎𝟗 + 𝟎. 𝟎𝟐 [ (𝟎. 𝟔𝟓𝟗𝒙𝟏𝟎. 𝟏𝟏𝟗𝟏𝟎𝟗𝒙 (𝟏 − ) + 𝒇(𝟎. 𝟎𝟑𝟑𝟑𝟑𝟑, 𝟏𝟎. 𝟏𝟗𝟖𝟗𝟑𝟒)]
𝟒 𝟗𝟗 𝟒
𝒀𝟏
𝒇(𝒙𝟏 , 𝒀𝟏 ) = 𝟎. 𝟔𝟓𝟗𝒀𝟏 (𝟏 − )
𝟗𝟗. 𝟎
𝟏 𝟑 𝟏𝟎. 𝟏𝟗𝟖𝟗𝟑𝟒
𝒀𝟏 = 𝟏𝟎 + 𝟎. 𝟎𝟐 [ (𝟓. 𝟗𝟖𝟔𝟖𝟖𝟒) + 𝒙(𝟎. 𝟔𝟓𝟗𝒙𝟏𝟎. 𝟏𝟗𝟖𝟗𝟑𝟒𝒙(𝟏 − )]
𝟒 𝟒 𝟗𝟗. 𝟎
𝟏 𝟑
𝒀𝟏 = 𝟏𝟎 + 𝟎. 𝟎𝟐 [ (𝟓. 𝟗𝟖𝟔𝟖𝟖𝟒) + 𝒙(𝟔. 𝟎𝟐𝟖𝟔𝟗𝟑)]
𝟒 𝟒
𝒀𝟏 = 𝟏𝟎. 𝟐𝟑𝟗𝟒𝟕𝟒
A modo de resumen se presenta la siguiente tabla:
𝒏 𝒙𝒏 𝒀𝒏
0 0.00 10.000000
1 0.02 10.119109
2 0.04 10.239474
18
Universidad Nacional de Salta
Facultad de Ingeniería
Ingeniería Química
Análisis Numérico
Verificación de la solución del problema con el programa en C++
Alumno: ARAMAYO, Juan Esteban – L.U: 314197 – Comisión: 2
19
