Iniciar sesiónCerrar
Correo:
Contraseña:
Entrar
Recordar sesión en este navegador
Iniciar sesiónIniciar sesiónCrear cuentaCrear cuenta
LWP » Foros » Matlab » Programa de ecuacion diferencial

Matlab - Programa de ecuacion diferencial

Volver
Nuevo Tema
<<>>
 
Vista:
sin imagen de perfil

Programa de ecuacion diferencial

Publicado por Patricio (6 intervenciones) el 01/07/2015 00:28:10
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
clear all,clf,clc
 
 
%Constantes

R=8.314; %ctte de los fases

To= 298; %temperatura inicial Kelvin

h=5; %velocidad de calentamiento Kelvin/s

 
%Constantes de velocidad para la descomposicion de biomasa

K1=(4.38E9)*exp(152700/(R*(To+(h*t))));%Descomposicion de la biomasa en gas

K2=(1.08E10)*exp(148000/(R*(To+(h*t))));%Descomposicion de la biomasa en bio-oil

K3=(3.375E10)*exp(111700/(R*(To+(h*t))));%Descomposicion de la biomasa en solido-intermediario

K4=(1.38E10)*exp(161000/(R*(To+(h*t))));%Descomposicion de solido intermediario en carbonizado

K5=(1.00E5)*exp(108000/(R*(To+(h*t))));%Descomposicion de bio-oil en carbonizado

K6=(4.28E6)*exp(108000/(R*(To+(h*t))));%Descomposicion de bio-oil en gas

 
%Programa de ecuacion diferencial

Eq=inline('[-(K1+K2+K3)*x(1);K3*x(1)-K4*x(2);K4*x(2)+K5*x(5);K1*x(1)+K6*x(5);K2*x(1)-(K6+K5)*x(5)]','t','x');
[t,x]=ode45(Eq,[0:1:70],[1 0 0 0 0]);

Las ecuaciones diferenciales son 5 y se pueden ver en Eq=inline, muchas gracias por su ayuda
Valora esta pregunta
Me gusta: Está pregunta es útil y esta claraNo me gusta: Está pregunta no esta clara o no es útil
0
Responder
Imágen de perfil de Daniel
Val: 269
Ha mantenido su posición en Matlab (en relación al último mes)
Gráfica de Matlab

Programa de ecuacion diferencial

Publicado por Daniel (264 intervenciones) el 01/07/2015 19:13:41
Supongo que tu pregunta es ¿porqué no anda?/¿como hago que funcione?

hay varios errores:
- las "constantes de velocidad ..." no son constantes sino que dependen de t, luego hay que definirlas como funciones, por ejemplo:

1
K1=@(t)(4.380E09)*exp(152700/(R*(To+(h*t))));


- la forma de definir el intervalo de resolución [0:1:70] no es válida, acá tenes dos opciones:
0:1:70 (sin corchetes) resuelve sólo en los puntos especificados
[0 70] resuelve en puntos que elije automaticamente (la cantidad que elije lo hace para garantizar una solucion adecuada)

- chequea que los valores de los K sean los correctos, especialmente el signo de la exponencial y el signo de la constante multiplicatiba, pues, si bien el código es correcto el ode no puede encontrar la solución. Aclaro que no estoy familiarizado que el proceso que estas ecuaciones están describiendo, pero creo que algo no esta bien.



acá te dejo el código con estas modificaciones



1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
clear all
clf,clc
 
%Constantes

R=8.314; %ctte de los fases

To= 298; %temperatura inicial Kelvin

h=5; %velocidad de calentamiento Kelvin/s

 
%Constantes de velocidad para la descomposicion de biomasa

K1=@(t)(4.380E09)*exp(152700/(R*(To+(h*t))));%Descomposicion de la biomasa en gas

K2=@(t)(1.080E10)*exp(148000/(R*(To+(h*t))));%Descomposicion de la biomasa en bio-oil

K3=@(t)(3.375E10)*exp(111700/(R*(To+(h*t))));%Descomposicion de la biomasa en solido-intermediario

K4=@(t)(1.380E10)*exp(161000/(R*(To+(h*t))));%Descomposicion de solido intermediario en carbonizado

K5=@(t)(1.000E05)*exp(108000/(R*(To+(h*t))));%Descomposicion de bio-oil en carbonizado

K6=@(t)(4.280E06)*exp(108000/(R*(To+(h*t))));%Descomposicion de bio-oil en gas

 
%Programa de ecuacion diferencial

Eq=@(t,x)([-(K1(t)+K2(t)+K3(t))*x(1);...
    K3(t)*x(1)-K4(t)*x(2);...
    K4(t)*x(2)+K5(t)*x(5);...
    K1(t)*x(1)+K6(t)*x(5);...
    K2(t)*x(1)-(K6(t)+K5(t))*x(5)]);
[t,x]=ode45(Eq,[0 70],[1 0 0 0 0]);
 
plot(t,x)


Saludos


Avisame como te fue

Daniel
Valora esta respuesta
Me gusta: Está respuesta es útil y esta claraNo me gusta: Está respuesta no esta clara o no es útil
1
Comentar
sin imagen de perfil

Programa de ecuacion diferencial

Publicado por Patricio (6 intervenciones) el 01/07/2015 23:05:18
Hola Daniel muchas gracias por tu ayuda la verdad me sirvió de bastante ya que soy Nuevo en esto de Matlab te cuento que trate de que el programa corra y la verdad que no encontró ninguna solución, revise los datos y serían los correctos. Ahora lo que hice es tratar de simplificar el sistema con un cambio de variable para usar solamente Temperatura.
De la sgte forma
T=f(t)
T=To+h*t ……(como en el primer programa que mande)
dT=h*dt………….. diferenciando
donde dt=dT/h

Originándose el siguiente programa siguiendo tus sugerencias
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
clear all
clf,clc
 
%Constantes

R=8.314; %ctte de los fases

h=5; %velocidad de calentamienT Kelvin/s

 
%Constantes de velocidad para la descomposicion de biomasa

K1=@(T)(4.380E09)*exp(152700/(R*T));%Descomposicion de la biomasa en gas

K2=@(T)(1.080E10)*exp(148000/(R*T));%Descomposicion de la biomasa en bio-oil

K3=@(T)(3.375E10)*exp(111700/(R*T));%Descomposicion de la biomasa en solido-intermediario

K4=@(T)(1.380E10)*exp(161000/(R*T));%Descomposicion de solido inTermediario en carbonizado

K5=@(T)(1.000E05)*exp(108000/(R*T));%Descomposicion de bio-oil en carbonizado

K6=@(T)(4.280E06)*exp(108000/(R*T));%Descomposicion de bio-oil en gas

 
%Programa de ecuacion diferencial

Eq=@(T,x)([(-(K1(T)+K2(T)+K3(T))*x(1))*1/h;...
    (K3(T)*x(1)-K4(T)*x(2))*1/h;...
    (K4(T)*x(2)+K5(T)*x(5))*1/h;...
    (K1(T)*x(1)+K6(T)*x(5))*1/h;...
    (K2(T)*x(1)-(K6(T)+K5(T))*x(5))*1/h]);
 
[T,x]=ode45(Eq,[298 900],[1 0 0 0 0]);
 
plot(T,x)

El problema ahora es que el la gráfica no grafica nada solo muestra una en blanco muchas gracias por tu ayuda, y gracias nuevamente de antemano y gracias por tu tiempo
Valora esta respuesta
Me gusta: Está respuesta es útil y esta claraNo me gusta: Está respuesta no esta clara o no es útil
0
Comentar
Imágen de perfil de Daniel
Val: 269
Ha mantenido su posición en Matlab (en relación al último mes)
Gráfica de Matlab

Programa de ecuacion diferencial

Publicado por Daniel (264 intervenciones) el 02/07/2015 02:06:51
te comento un poco mas de la lógica que se usa para resolver estas ecuaciones en forma numerica para este sistema particular:

la idea básica detrás de resolver una ecuación diferencial en forma numérica es la siguiente:

supongamos para no dar tantas vueltas que tenemos solo la primer ecuación del sistema:

x'=-(k1+k2+k3)*x

o bien dx/dT=-(k1+k2+k3)*x

si tomo un DT es suficientemente pequeño, dDx/dT será aproximadamente (xn-xv)/DT, donde xn es el valor en t+DT y xv es el valor en t, por lo cual xn-xv=-(k1+k2+k3)*xv*DT (donde el x lo escribí como xv ya que la diferencia entre xv y xn quiero que sea pequeña).

entonces xn=xv-(k1+k2+k3)*xv*DT o bien xn=xv(1-(k1+k2+k3)*DT)

como dije antes quiero que la solucion vaya variando de a poco de manera que xv y xn sean muy parecidos, luego esto impone la condición (k1+k2+k3)*DT<<1.

en tu caso k1+k2+k3 es del orden de 1e36, con lo cual el DT tiene que ser mucho menor a 1e-36, al menos al principio (luego los k van variando).

por eso al matlab le resulta imposible resolver este problema para intervalos de temperatura tan grande. Como yo creo que no estas queriendo describir un sistema "físico" que varíe tan fuertemente con la temperatura es que te decía que revises los números y unidades de los k's.

yo probé con intervalos de temperatura muy pequeños [0 4e-36] y obtuve resultados

1
[t,x]=ode45(Eq,[0 1e-34],[1 0 0 0 0]);

Te dejo que lo veas un poco mas de tu parte y cualquier cosa seguimos en contacto por aquí, si te puedo seguir dando una mano.

Saludos

Daniel
Valora esta respuesta
Me gusta: Está respuesta es útil y esta claraNo me gusta: Está respuesta no esta clara o no es útil
0
Comentar