PROYECTO DE ESTUDIO DE CASOS: ECUACIONES DIFERENCIALES
PARCIALES
CAPITULO 5
EJERCICIO 5.1:
Flujo dinámico en un tanque puede ser modelado haciendo un balance de masa sobre el
fluido en el tanque (Figura 5.10). La naturaleza del resultado ODE-IVP depende del
modelo utilizado para la válvula de salida. Si una válvula lineal es
A = área sección transversal
FIGURA 5.10: Flujo dinámico en un esquema de tanque
Asumido, entonces el modelo ODE es lineal. Una descripción más precisa de la válvula
hace que el flujo de salida sea una función no lineal de la altura. Ambos casos se
consideran en lo que sigue.
Balance de masa:
𝑑(𝑚𝑎𝑠𝑎 𝑒𝑛 𝑡𝑎𝑛𝑞𝑢𝑒 𝑎 𝑐𝑢𝑎𝑙𝑞𝑢𝑖𝑒𝑟 𝑡𝑖𝑒𝑚𝑝𝑜)
= 𝑒𝑛𝑡𝑟𝑎𝑑𝑎 − 𝑠𝑎𝑙𝑖𝑑𝑎 + 𝑔𝑒𝑛𝑒𝑟𝑎𝑐𝑖𝑜𝑛
𝑑𝑡
Asuma un componente puro (único) de fluido de densidad constante:
𝑑𝜌𝑉
= 𝐹𝑒𝑛𝑡𝑟𝑎 − 𝐹𝑠𝑎𝑙𝑒
𝑑𝑡
Donde:
𝝆 = 𝒅𝒆𝒏𝒔𝒊𝒅𝒂𝒅
𝑽 = 𝒗𝒐𝒍𝒖𝒎𝒆𝒏
𝑭 = 𝒓𝒂𝒛𝒐𝒏 𝒅𝒆 𝒇𝒍𝒖𝒋𝒐 𝒅𝒆 𝒎𝒂𝒔𝒂
En términos de área transversal y altura,
𝑑𝜌𝐴ℎ
= 𝐹𝑒𝑛𝑡𝑟𝑎 − 𝐹𝑠𝑎𝑙𝑒
𝑑𝑡
donde la densidad y el área son constantes.
𝑑ℎ 𝐹𝑒𝑛𝑡𝑟𝑎 − 𝐹𝑠𝑎𝑙𝑒
=
𝑑𝑡 𝜌𝐴
𝑭𝒆𝒏𝒕𝒓𝒂 puede ser:
a. Directamente proporcional a h: 𝑭𝒔𝒂𝒍𝒆 = 𝑪, 𝒉 en cuyo caso la ecuación del modelo
es:
𝑑ℎ 𝐹𝑒𝑛𝑡𝑟𝑎 − 𝐶, ℎ
=
𝑑𝑡 𝜌𝐴
𝟏⁄
b. Una función de h, como 𝑭𝒔𝒂𝒍𝒆 = 𝑪, 𝒉 𝟐 y el modelo se convierte en:
𝑑ℎ 𝐹𝑒𝑛𝑡𝑟𝑎 − 𝐶, √ℎ
=
𝑑𝑡 𝜌𝐴
Los datos adicionales pertinentes son los siguientes:
Diámetro del tanque = 𝟔 𝒇𝒕
Densidad del líquido = 𝟔𝟐. 𝟓 𝒍𝒃/𝒇𝒕𝟑
Constante de la válvula 𝑪𝒗 = 𝟐𝟓𝟎 𝒍𝒃/𝒉 − 𝒇𝒕 𝒑𝒂𝒓𝒂 𝒆𝒍 𝒄𝒂𝒔𝒐 𝒍𝒊𝒏𝒆𝒂𝒍.
Constante de la Válvula 𝑪𝒗 = 𝟓𝟎𝟎 𝒍𝒃/𝒉 − 𝒇𝒕𝟎.𝟓 𝒑𝒂𝒓𝒂 𝒆𝒍 𝒄𝒂𝒔𝒐 𝒏𝒐 𝒍𝒊𝒏𝒆𝒂𝒍.
Altura inicial = 𝟒 𝒇𝒕
𝑭𝒆𝒏𝒕𝒓𝒂 = 𝟏𝟎𝟎𝟎 𝒍𝒃/𝒉 (Flujo inicial de entrada en estado estacionario)
𝒕 = 𝒕𝒊𝒆𝒎𝒑𝒐(𝒉)
𝑭𝒆𝒏𝒕𝒓𝒂 sufre un cambio en el tiempo 0+ a 𝟏𝟒𝟎𝟎 𝒍𝒃/𝒉 hasta el tiempo 10h, después de lo
cual retorna a 𝟏𝟎𝟎𝟎 𝒍𝒃/𝒉.
a. Resolver el lineal ODE (𝑭𝒔𝒂𝒍𝒆 𝒑𝒓𝒐𝒑𝒐𝒓𝒄𝒊𝒐𝒏𝒂𝒍 𝒂 𝒉) usando el método de Euler,
empezar con una ∆𝒕 𝒅𝒆 𝟎. 𝟓 𝒚 𝒅𝒆𝒏𝒖𝒆𝒗𝒐 𝒄𝒐𝒏 ∆𝒕 = 𝟎. 𝟎𝟓 para determinar una ∆𝒕.
Use IF function para generar los valores adecuados para 𝑭𝒆𝒏𝒕𝒓𝒂.
𝟏⁄
b. Repita la parte a para el caso no lineal (𝑭𝒔𝒂𝒍𝒆 𝒑𝒓𝒐𝒑𝒐𝒓𝒄𝒊𝒐𝒏𝒂𝒍 𝒂 𝒉 𝟐 ).
Figura del problema:
Análisis y ecuaciones:
𝐹𝑒𝑛𝑡𝑟𝑎 puede ser:
Directamente proporcional a h: 𝐹𝑠𝑎𝑙𝑒 = 𝐶, ℎ en cuyo caso la ecuación del modelo es:
𝑑ℎ 𝐹𝑒𝑛𝑡𝑟𝑎 − 𝐶, ℎ
=
𝑑𝑡 𝜌𝐴
1⁄
Una función de h, como 𝐹𝑠𝑎𝑙𝑒 = 𝐶, ℎ 2 y el modelo se convierte en:
𝑑ℎ 𝐹𝑒𝑛𝑡𝑟𝑎 − 𝐶, √ℎ
=
𝑑𝑡 𝜌𝐴
Para (a):
Código Matlab:
% Ejercicio 5.1 libro chNumerical Methods for Chemical Engineer.
clear all
clc
disp('**************************************************************')
disp('* UNIVERSIDAD NACIONAL DE TRUJILLO *')
disp('* FACULTAD DE INGENIERIA QUIMICA *')
disp('* ESCUELA PROFESIONAL DE INGENIERA QUIMICA *')
disp('*CURSO: METODOS NUMERICOS *')
disp('*Libro: Numerical Methods for Chemical Engineer.- capitulo 5 *')
disp('*DOCENTE: Dr. Guillermo Evangelista Benites *')
disp('*Aumno: Javier Loredo Luna Victoria *')
disp('**************************************************************')
disp(' ')
%% INGRESO DE DATOS1
xf=1; % valor al cual se desea conocer la func.
h=0.5;
xo=0; yo=4; %condicion inicial
i=0;
fprintf('\n\n\n\t\t\t\t MÉTODO DE RUNGE-KUTTA (RK-4)')
fprintf('\n\t\t\t\t===================================\n')
fprintf('\n\t\t\t\t\tpara dt = 0.5\n')
fprintf('\t\t\t\t------------------------------\n')
fprintf('\t\t\t \ti \t t(i) \tH(i) \n')
fprintf('\t\t\t\t------------------------------\n')
fprintf('\t\t\t \t%2d \t%6.4f \t%8.6f\n',i,xo,yo)
N=(xf-xo)/h;
if xf==0
f=inline('0*x+1000/(3.1416*62.5*9)-
250/(3.1416*62.5*9)*y','x','y');
elseif xf>0 & xf<=10
f=inline('0*x+1400/(3.1416*62.5*9)-
250/(3.1416*62.5*9)*y','x','y');
elseif xf>10
f=inline('0*x+1000/(3.1416*62.5*9)-
250/(3.1416*62.5*9)*y','x','y');
end
x1=zeros(1,length(N));
x1(1,1)=xo;
y1=zeros(1,length(N));
y1(1,1)=yo;
while i<N
k1=f(xo,yo);
k2=f(xo+h/2,yo+h*k1/2);
k3=f(xo+h/2,yo+h*k2/2);
k4=f(xo+h,yo+h*k3);
yo=yo+h/6*(k1+2*k2+2*k3+k4);
xo=xo+h;
i=i+1;
fprintf('\t\t\t \t%2d \t%6.4f \t%8.6f\n',i,xo,yo)
x1(1,i+1)=xo;
y1(1,i+1)=yo;
end
fprintf('\t\t\t\t-----------------------------\n')
xx=0:h:xf;
[T,Y]=ode45(f,xx,4);
x_analit=Y(end);
ERP=(abs(yo-x_analit))/x_analit*100;
fprintf('\t\t\t\t El error = %6.4e\n',ERP)
subplot(2,2,1);
plot(x1,y1,'ro-')
xlabel('t'); ylabel('H')
title('dt=0.5');grid
hold on
plot(T,Y,'b')
legend('Runge-kuta','ode45')
%=======================================
xf=1; % valor al cual se desea conocer la func.
h=0.05;
xo=0; yo=4; %condicion inicial
i=0;
fprintf('\n\n\t\t\t\t\tpara dt = 0.05\n')
fprintf('\t\t\t\t------------------------------\n')
fprintf('\t\t\t \ti \t t(i) \tH(i) \n')
fprintf('\t\t\t\t------------------------------\n')
fprintf('\t\t\t \t%2d \t%6.4f \t%8.6f\n',i,xo,yo)
N=(xf-xo)/h;
if xf==0
f=inline('0*x+1000/(3.1416*62.5*9)-
250/(3.1416*62.5*9)*y','x','y');
elseif xf>0 & xf<=10
f=inline('0*x+1400/(3.1416*62.5*9)-
250/(3.1416*62.5*9)*y','x','y');
elseif xf>10
f=inline('0*x+1000/(3.1416*62.5*9)-
250/(3.1416*62.5*9)*y','x','y');
end
x1=zeros(1,length(N));
x1(1,1)=xo;
y1=zeros(1,length(N));
y1(1,1)=yo;
while i<N
k1=f(xo,yo);
k2=f(xo+h/2,yo+h*k1/2);
k3=f(xo+h/2,yo+h*k2/2);
k4=f(xo+h,yo+h*k3);
yo=yo+h/6*(k1+2*k2+2*k3+k4);
xo=xo+h;
i=i+1;
fprintf('\t\t\t \t%2d \t%6.4f \t%8.6f\n',i,xo,yo)
x1(1,i+1)=xo;
y1(1,i+1)=yo;
end
fprintf('\t\t\t\t-----------------------------\n')
xx=0:h:xf;
[T,Y]=ode45(f,xx,4);
x_analit=Y(end);
ERP=(abs(yo-x_analit))/x_analit*100;
fprintf('\t\t\t\t El error = %6.4e\n',ERP)
subplot(2,2,2);
plot(x1,y1,'ro-')
xlabel('t'); ylabel('H')
title('dt=0.05');grid
hold on
plot(T,Y,'b')
legend('Runge-kuta','ode45')
%=================================================
xf=1; % valor al cual se desea conocer la func.
h=0.05;
xo=0; yo=4; %condicion inicial
i=0;
fprintf('\n\n\t\t\t\t\tpara dt = 0.005\n')
fprintf('\t\t\t\t------------------------------\n')
fprintf('\t\t\t \ti \t t(i) \tH(i) \n')
fprintf('\t\t\t\t------------------------------\n')
fprintf('\t\t\t \t%2d \t%6.4f \t%8.6f\n',i,xo,yo)
N=(xf-xo)/h;
if xf==0
f=inline('0*x+1000/(3.1416*62.5*9)-
250/(3.1416*62.5*9)*y','x','y');
elseif xf>0 & xf<=10
f=inline('0*x+1400/(3.1416*62.5*9)-
250/(3.1416*62.5*9)*y','x','y');
elseif xf>10
f=inline('0*x+1000/(3.1416*62.5*9)-
250/(3.1416*62.5*9)*y','x','y');
end
x1=zeros(1,length(N));
x1(1,1)=xo;
y1=zeros(1,length(N));
y1(1,1)=yo;
while i<N
k1=f(xo,yo);
k2=f(xo+h/2,yo+h*k1/2);
k3=f(xo+h/2,yo+h*k2/2);
k4=f(xo+h,yo+h*k3);
yo=yo+h/6*(k1+2*k2+2*k3+k4);
xo=xo+h;
i=i+1;
fprintf('\t\t\t \t%2d \t%6.4f \t%8.6f\n',i,xo,yo)
x1(1,i+1)=xo;
y1(1,i+1)=yo;
end
fprintf('\t\t\t\t-----------------------------\n')
xx=0:h:xf;
[T,Y]=ode45(f,xx,4);
x_analit=Y(end);
ERP=(abs(yo-x_analit))/x_analit*100;
fprintf('\t\t\t\t El error = %6.4e\n',ERP)
subplot(2,2,3);
plot(x1,y1,'ro-')
xlabel('t'); ylabel('H')
title('dt=0.005');grid
hold on
plot(T,Y,'b')
legend('Runge-kuta','ode45')
disp('El valor para dt = 0.005')
Resultados:
Usando el software MATLAB para resolver sistemas de ecuaciones diferenciales, obtenemos:
Graficas
Para (b):
Código Polymath:
A = 3.1416 * D ^ 2 / 4 # area
D = 6 # diametro del cilindro
t(0) = 0
t(f) = 0.85
Cv = 250 # valor para enunciado a
F = 1400 # para el valor de 0<t<=10
p = 62.5 # densidad del fluido
d(h)/d(t) = (F - Cv * sqrt(h)) / (p * A) # ecuacion para enunciado a
h(0) = 4
Resultados:
Usando el software Polymath para resolver sistemas de ecuaciones diferenciales se obtienen los
siguientes resultados:
Código Matlab:
% Ejercicio 5.1 libro chNumerical Methods for Chemical Engineer.
clear all
clc
disp('**************************************************************')
disp('* UNIVERSIDAD NACIONAL DE TRUJILLO *')
disp('* FACULTAD DE INGENIERIA QUIMICA *')
disp('* ESCUELA PROFESIONAL DE INGENIERA QUIMICA *')
disp('*CURSO: METODOS NUMERICOS *')
disp('*Libro: Numerical Methods for Chemical Engineer.- capitulo 5 *')
disp('*DOCENTE: Dr. Guillermo Evangelista Benites *')
disp('*Aumno: Javier Loredo Luna Victoria *')
disp('**************************************************************')
disp(' ')
%% INGRESO DE DATOS1
xf=1; % valor al cual se desea conocer la func.
h=0.5;
xo=0; yo=4; %condicion inicial
i=0;
fprintf('\n\n\n\t\t\t\t MÉTODO DE RUNGE-KUTTA (RK-4)')
fprintf('\n\t\t\t\t===================================\n')
fprintf('\t\t\t\t\tpara dt = 0.5\n')
fprintf('\n\t\t\t\t------------------------------\n')
fprintf('\t\t\t \ti \t t(i) \tH(i) \n')
fprintf('\t\t\t\t------------------------------\n')
fprintf('\t\t\t \t%2d \t%6.4f \t%8.6f\n',i,xo,yo)
N=(xf-xo)/h;
if xf==0
f=inline('0*x+1000/(3.1416*62.5*9)-
500/(3.1416*62.5*9)*y^0.5','x','y');
elseif xf>0 & xf<=10
f=inline('0*x+1400/(3.1416*62.5*9)-
500/(3.1416*62.5*9)*y^0.5','x','y');
elseif xf>10
f=inline('0*x+1000/(3.1416*62.5*9)-
500/(3.1416*62.5*9)*y^0.5','x','y');
end
x1=zeros(1,length(N));
x1(1,1)=xo;
y1=zeros(1,length(N));
y1(1,1)=yo;
while i<N
k1=f(xo,yo);
k2=f(xo+h/2,yo+h*k1/2);
k3=f(xo+h/2,yo+h*k2/2);
k4=f(xo+h,yo+h*k3);
yo=yo+h/6*(k1+2*k2+2*k3+k4);
xo=xo+h;
i=i+1;
fprintf('\t\t\t \t%2d \t%6.4f \t%8.6f\n',i,xo,yo)
x1(1,i+1)=xo;
y1(1,i+1)=yo;
end
fprintf('\t\t\t\t-----------------------------\n')
xx=0:h:xf;
[T,Y]=ode45(f,xx,4);
x_analit=Y(end);
ERP=(abs(yo-x_analit))/x_analit*100;
fprintf('\t\t\t\t El error = %6.4e\n',ERP)
subplot(1,2,1);
plot(x1,y1,'ro-')
xlabel('t'); ylabel('H')
title('dt=0.5');grid
hold on
plot(T,Y,'b')
legend('Runge-kuta','ode45')
%=======================================
xf=1; % valor al cual se desea conocer la func.
h=0.05;
xo=0; yo=4; %condicion inicial
i=0;
fprintf('\n\t\t\t\t\tpara dt = 0.05\n')
fprintf('\t\t\t\t------------------------------\n')
fprintf('\t\t\t \ti \t t(i) \tH(i) \n')
fprintf('\t\t\t\t------------------------------\n')
fprintf('\t\t\t \t%2d \t%6.4f \t%8.6f\n',i,xo,yo)
N=(xf-xo)/h;
if xf==0
f=inline('0*x+1000/(3.1416*62.5*9)-
500/(3.1416*62.5*9)*y^0.5','x','y');
elseif xf>0 & xf<=10
f=inline('0*x+1400/(3.1416*62.5*9)-
500/(3.1416*62.5*9)*y^0.5','x','y');
elseif xf>10
f=inline('0*x+1000/(3.1416*62.5*9)-
500/(3.1416*62.5*9)*y^0.5','x','y');
end
x1=zeros(1,length(N));
x1(1,1)=xo;
y1=zeros(1,length(N));
y1(1,1)=yo;
while i<N
k1=f(xo,yo);
k2=f(xo+h/2,yo+h*k1/2);
k3=f(xo+h/2,yo+h*k2/2);
k4=f(xo+h,yo+h*k3);
yo=yo+h/6*(k1+2*k2+2*k3+k4);
xo=xo+h;
i=i+1;
fprintf('\t\t\t \t%2d \t%6.4f \t%8.6f\n',i,xo,yo)
x1(1,i+1)=xo;
y1(1,i+1)=yo;
end
fprintf('\t\t\t\t-----------------------------\n')
xx=0:h:xf;
[T,Y]=ode45(f,xx,4);
x_analit=Y(end);
ERP=(abs(yo-x_analit))/x_analit*100;
fprintf('\t\t\t\t El error = %6.4e\n',ERP)
subplot(1,2,2);
plot(x1,y1,'ro-')
xlabel('t'); ylabel('H')
title('dt=0.05');grid
hold on
plot(T,Y,'b')
legend('Runge-kuta','ode45')
disp('el valor adecuado para dt = 0.05')
Resultados:
Usando el software MATLAB para resolver sistemas de ecuaciones diferenciales se obtienen los
siguientes resultados:
Graficas
EJERCICIO 5.2:
Resolver el ODEs del ejercicio 5.1 usando el método de Runge-Kutta de segundo orden.
Encontrar valores aceptables para h cuando se usa el método de Runge-Kutta de
segundo orden.
Solución
Para hallar el valor correcto de h se probó valores de 0.03, 0.02, 0.01, para luego
comparar con ODE23 de MatLab. Dando como valor correcto h=0.02 ya que dio la solución
del problema en menos iteraciones, dado que el algoritmo de Runge-Kutta de segundo
orden genera más iteraciones que el de cuarto orden para poder obtener el valor de la
solución de las EDOs.
%% INGRESO DE DATOS1
%Prueba con h=0.03
xf=1; % valor al cual se desea conocer la func.
h=0.03;%Variando h
xo=0; yo=4; %condicion inicial
i=0;
fprintf('\n\n\n\t\t\t\t MÉTODO DE RUNGE-KUTTA (RK-2)')
fprintf('\n\t\t\t\t===================================\n')
fprintf('\n\t\t\t\t\tpara dt = 0.03\n')
fprintf('\t\t\t\t------------------------------\n')
fprintf('\t\t\t \ti \t t(i) \tH(i) \n')
fprintf('\t\t\t\t------------------------------\n')
fprintf('\t\t\t \t%2d \t%6.4f \t%8.6f\n',i,xo,yo)
N=(xf-xo)/h;
%Reemplazando los valores para obtener las funciones en inline()
if xf==0
f=inline('0*x+1000/(3.1416*62.5*9)-
250/(3.1416*62.5*9)*y','x','y');
elseif xf>0 & xf<=10
f=inline('0*x+1400/(3.1416*62.5*9)-
250/(3.1416*62.5*9)*y','x','y');
elseif xf>10
f=inline('0*x+1000/(3.1416*62.5*9)-
250/(3.1416*62.5*9)*y','x','y');
end
x1=zeros(1,length(N));
x1(1,1)=xo;
y1=zeros(1,length(N));
y1(1,1)=yo;
while i<N
k1=f(xo,yo);
k2=f(xo+h,yo+k1*h);
yo=yo+h/2*(k1+k2);
xo=xo+h;
i=i+1;
fprintf('\t\t\t \t%2d \t%6.4f \t%8.6f\n',i,xo,yo)
x1(1,i+1)=xo;
y1(1,i+1)=yo;
end
fprintf('\t\t\t\t-----------------------------\n')
xx=0:h:xf;
[T,Y]=ode23(f,xx,4);
x_analit=Y(end);
ERP=(abs(yo-x_analit))/x_analit*100;
fprintf('\t\t\t\t El error = %6.4e\n',ERP)
subplot(2,2,1);
plot(x1,y1,'ro-')
xlabel('t'); ylabel('H')
title('dt=0.03');grid
hold on
plot(T,Y,'b')
legend('Runge-kuta','ode23')
%=======================================
%Prueba con h=0.02
xf=1; % valor al cual se desea conocer la func.
h=0.02;%Variando h
xo=0; yo=4; %condicion inicial
i=0;
fprintf('\n\n\t\t\t\t\tpara dt = 0.02\n')
fprintf('\t\t\t\t------------------------------\n')
fprintf('\t\t\t \ti \t t(i) \tH(i) \n')
fprintf('\t\t\t\t------------------------------\n')
fprintf('\t\t\t \t%2d \t%6.4f \t%8.6f\n',i,xo,yo)
N=(xf-xo)/h;
if xf==0
f=inline('0*x+1000/(3.1416*62.5*9)-
250/(3.1416*62.5*9)*y','x','y');
elseif xf>0 & xf<=10
f=inline('0*x+1400/(3.1416*62.5*9)-
250/(3.1416*62.5*9)*y','x','y');
elseif xf>10
f=inline('0*x+1000/(3.1416*62.5*9)-
250/(3.1416*62.5*9)*y','x','y');
end
x1=zeros(1,length(N));
x1(1,1)=xo;
y1=zeros(1,length(N));
y1(1,1)=yo;
while i<N
k1=f(xo,yo);
k2=f(xo+h,yo+k1*h);
yo=yo+h/2*(k1+k2);
xo=xo+h;
i=i+1;
fprintf('\t\t\t \t%2d \t%6.4f \t%8.6f\n',i,xo,yo)
x1(1,i+1)=xo;
y1(1,i+1)=yo;
end
fprintf('\t\t\t\t-----------------------------\n')
xx=0:h:xf;
[T,Y]=ode23(f,xx,4);
x_analit=Y(end);
ERP=(abs(yo-x_analit))/x_analit*100;
fprintf('\t\t\t\t El error = %6.4e\n',ERP)
subplot(2,2,2);
plot(x1,y1,'ro-')
xlabel('t'); ylabel('H')
title('dt=0.02');grid
hold on
plot(T,Y,'b')
legend('Runge-kuta','ode23')
%=================================================
%Prueba con h=0.01
xf=1; % valor al cual se desea conocer la func.
h=0.01; %Variando h
xo=0; yo=4; %condicion inicial
i=0;
fprintf('\n\n\t\t\t\t\tpara dt = 0.01\n')
fprintf('\t\t\t\t------------------------------\n')
fprintf('\t\t\t \ti \t t(i) \tH(i) \n')
fprintf('\t\t\t\t------------------------------\n')
fprintf('\t\t\t \t%2d \t%6.4f \t%8.6f\n',i,xo,yo)
N=(xf-xo)/h;
if xf==0
f=inline('0*x+1000/(3.1416*62.5*9)-
250/(3.1416*62.5*9)*y','x','y');
elseif xf>0 & xf<=10
f=inline('0*x+1400/(3.1416*62.5*9)-
250/(3.1416*62.5*9)*y','x','y');
elseif xf>10
f=inline('0*x+1000/(3.1416*62.5*9)-
250/(3.1416*62.5*9)*y','x','y');
end
x1=zeros(1,length(N));
x1(1,1)=xo;
y1=zeros(1,length(N));
y1(1,1)=yo;
while i<N
k1=f(xo,yo);
k2=f(xo+h,yo+k1*h);
yo=yo+h/2*(k1+k2);
xo=xo+h;
i=i+1;
fprintf('\t\t\t \t%2d \t%6.4f \t%8.6f\n',i,xo,yo)
x1(1,i+1)=xo;
y1(1,i+1)=yo;
end
fprintf('\t\t\t\t-----------------------------\n')
xx=0:h:xf;
[T,Y]=ode23(f,xx,4);
x_analit=Y(end);
ERP=(abs(yo-x_analit))/x_analit*100;
fprintf('\t\t\t\t El error = %6.4e\n',ERP)
subplot(2,2,3);
plot(x1,y1,'ro-')
xlabel('t'); ylabel('H')
title('dt=0.01');grid
hold on
plot(T,Y,'b')
legend('Runge-kuta','ode23')
disp('El valor para dt = 0.005')
Resultados del Algoritmo:
*Con h=0.03 =dt
# iteraciones= 34
Valor Final: 4.214995
Error: 1.3993e-01
*Con h=0.03 = dt
# iteraciones= 50
Valor Final: 4.211071
Error: 6.7630e-06
*Con h=0.03 =dt
# iteraciones= 100
Valor Final: 4.211071
Error: 2.0641e-06
EJERCICIO 5.3:
Resolver el ODEs del ejercicio 5.1 usando el programa de Euler VBA como dan en el
ejempo 5.7.
fprintf('\n\n\n\t\t\t\t MÉTODO DE EULER ')
clear all; clc
syms x
syms y
x=0;xf=1;y=4;
h=0.05;
n=(xf-x)/h;
if xf==0
f=inline('0*x+1000/(3.1416*62.5*9)-
250/(3.1416*62.5*9)*y','x','y');
elseif xf>0 & xf<=10
f=inline('0*x+1400/(3.1416*62.5*9)-
250/(3.1416*62.5*9)*y','x','y');
elseif xf>10
f=inline('0*x+1000/(3.1416*62.5*9)-
250/(3.1416*62.5*9)*y','x','y');
end
disp(' x y');
for i=1:n+1 %euler
y1=feval(f,x,y);
hy1=h*y1;
fprintf('\n%0.1f %0.4f ',x,y);
y=y+hy1;
x=x+h;
end
EJERCICIO 5.4:
Supongamos que las siguientes reacciones químicas tienen lugar en un reactor continuo
de tanque agitado (CSTR):
𝒌𝟏
𝑨⇔𝑩
Donde las constantes de velocidad son las siguientes:
𝒌𝟏 = 𝟏 𝒎𝒊𝒏−𝟏 , 𝒌𝟐 = 𝟎. 𝟓 𝒎𝒊𝒏−𝟏
La carga al reactor es todo el componente A, así que inicialmente, las concentraciones
dentro del reactor son
𝑪𝑨𝟎 = 𝟏 𝑪𝑩𝟎 = 𝟎 (𝒈𝒎𝒐𝒍/𝑳)
Un balance de masa de estado inestable en cada componente conduce al siguiente
conjunto de ODEs:
𝒅𝑪𝑨
= −𝒌𝟏 𝑪𝑨 + 𝒌𝟐 𝑪𝑩
𝒅𝒕
𝒅𝑪𝑩
= 𝒌𝟏 𝑪𝑨 − 𝒌𝟐 𝑪𝑩
𝒅𝒕
Resolver este sistema ODE-IVP utilizando el método de Euler en Excel con un paso de
tiempo de 0,1 min y repetir con un paso de tiempo de 0,01 min. Si se produce una
diferencia significativa (comparar los valores en un momento específico para ver si son
diferentes), repetir a con un paso de tiempo de 0,001, etc. hasta que se encuentre un
paso de tiempo satisfactorio. Integrar a tiempo = 2 min. Graficar la solución en cada
caso y compararlas. Reporte la diferencia porcentual en las soluciones a tiempo = 1 min.
Realice experimentos cambiando las constantes de velocidad en una de las soluciones
para ver cómo Excel recalcula inmediatamente las soluciones y actualiza la gráfica.
Figura del problema:
Análisis y ecuaciones:
Un balance de masa de estado inestable en cada componente conduce al siguiente
conjunto de ODEs:
𝒅𝑪𝑨
= −𝒌𝟏 𝑪𝑨 + 𝒌𝟐 𝑪𝑩
𝒅𝒕
𝒅𝑪𝑩
= 𝒌𝟏 𝑪𝑨 − 𝒌𝟐 𝑪𝑩
𝒅𝒕
Código Matlab:
function z = E54(f,g,CA,CB,t,h)
u1=h*eval('f(t,CA,CB)');
v1=h*eval('g(t,CA,CB)');
u2=h*eval('f(t+h/2,CA+u1/2,CB+v1/2)');
v2=h*eval('g(t+h/2,CA+u1/2,CB+v1/2)');
u3=h*eval('f(t+h/2,CA+u2/2,CB+v2/2)');
v3=h*eval('g(t+h/2,CA+u2/2,CB+v2/2)');
u4=h*eval('f(t+h,CA+u3,CB+v3)');
v4=h*eval('g(t+h,CA+u3,CB+v3)');
f1=CA+(1/6)*(u1+2*u2+2*u3+u4);
g1=CB+(1/6)*(v1+2*v2+2*v3+v4);
z=[f1,g1];
% Ejercicio 5.4 libro chNumerical Methods for Chemical Engineer.
clear all
clc
disp('**************************************************************')
disp('* UNIVERSIDAD NACIONAL DE TRUJILLO *')
disp('* FACULTAD DE INGENIERIA QUIMICA *')
disp('* ESCUELA PROFESIONAL DE INGENIERA QUIMICA *')
disp('*CURSO: METODOS NUMERICOS *')
disp('*Libro: Numerical Methods for Chemical Engineer.- capitulo 5 *')
disp('*DOCENTE: Dr. Guillermo Evangelista Benites *')
disp('*Aumno: Javier Loredo Luna Victoria *')
disp('**************************************************************')
disp(' ')
A=inline('-CA+(0.5*CB)','t','CA','CB');
B=inline('CA-(0.5*CB)','t','CA','CB');
t0=0;
tf=2;
h=0.01;
t=t0;
CA0=1.0;
CB0=0;
m = E54(A,B,CA0,CB0,t,h);
tab = [t,m];
n=m;
while(t<tf)
CA0=m(1);
CB0=m(2);
t=t+h;
tab=[tab;t,CA0,CB0];
m = E54(A,B,CA0,CB0,t,h);
n = [n;m];
end
tab;
fprintf('\t--------------------------------\n')
fprintf('\t\t t CA CB\n')
fprintf('\t--------------------------------\n')
fprintf('\t%8.2f %10.5f %10.5f\n',tab')
fprintf('\t--------------------------------\n')
u=0:h:tf;
plot(u',n)
title('CONCENTRACION vs TIEMPO','color','b')
xlabel('Tiempo, t','color','b')
ylabel('Concentracion, C,','color','b')
text(0.35,0.8,'CA','color','r');text(0.3,0.35,'CB','color','r
');
legend('CA: con paso de t= 0.01min','CB: con paso de t =
0.01min')
grid on
Resultados:
Usando el software MATLAB para resolver sistemas de ecuaciones diferenciales se
obtienen los siguientes resultados:
Grafica
Con paso de t = 0.01min.
𝐶𝐴 = 0.99007 𝑔𝑚𝑜𝑙/𝐿
𝐶𝐵 = 0.00993 𝑔𝑚𝑜𝑙/𝐿
Grafica
EJERCICIO 5.5:
Resuelva los EDO del ejercicio 5.4 utilizando el método de segundo orden Runge-Kutta.
Compare los resultados para ∆𝒕 = 0.1 min y para ∆𝒕 = 0.01 min.
Solución:
Matlab tiene incluido el método Runge-Kutta de segundo orden en la función Ode23,
entonces haremos uso de esa función:
*para ∆𝒕 = 0.01 min
clear all,clc
% k1=1;
% k2=0.5;
Eq=inline('[-1*x(1)+0.5*x(2);1*x(1)-0.5*x(2)]','t','x'); %ecuaciones de
diseño de reactor
[t,x]=ode23(Eq,[0:0.01:2],[1 0]); %dt=0.01
disp([t,x])
plot(t,x(:,1),'-r',t,x(:,2),'-b','LineWidth',2),grid on
title('Solución de un sistema de ecuaciones reactivas','Fontsize',15);
text(0.2,0.8,'CA');text(0.8,0.1,'CB');
xlabel('Tiempo');ylabel('Concentración');
*para ∆𝒕 = 0.1 min
clear all,clc
% k1=1;
% k2=0.5;
Eq=inline('[-1*x(1)+0.5*x(2);1*x(1)-0.5*x(2)]','t','x'); %ecuaciones de
diseño de reactor
[t,x]=ode23(Eq,[0:0.1:2],[1 0]); %dt=0.1
disp([t,x])
plot(t,x(:,1),'-r',t,x(:,2),'-b','LineWidth',2),grid on
title('Solución de un sistema de ecuaciones reactivas','Fontsize',15);
text(0.2,0.8,'CA');text(0.8,0.1,'CB');
xlabel('Tiempo');ylabel('Concentración');
Resultados:
para ∆𝒕 = 0.01 min
…..
para ∆𝒕 = 0.1 min
PARA t = 0.01min.
𝐶𝐴 = 0.9901 𝑔𝑚𝑜𝑙/𝐿
𝐶𝐵 = 0.0099 𝑔𝑚𝑜𝑙/𝐿
PARA t = 0.1min.
𝐶𝐴 = 0.9071 𝑔𝑚𝑜𝑙/𝐿
𝐶𝐵 = 0.0929 𝑔𝑚𝑜𝑙/𝐿
EJERCICIO 5.6:
Resuelva los EDO del ejercicio 5.4 utilizando el programa de Euler VBA
como dan en el ejemplo 5.7. Use los pasos de tiempo de 0.1, 0.01 y 0.001
minutos, respectivamente y deducir cuál es el más apropiado para una
solución aceptable.
Solución:
fprintf('\n\n\n\t\t\t\t MÉTODO DE EULER ')
clear all; clc
syms x
syms y
x=0;xf=1;y=1;
h=0.1;
n=(xf-x)/h;
f=inline('-1*x+0.5*y','x','y');
disp('CONCENTRACION CA')
disp(' x y');
for i=1:n+1
y1=feval(f,x,y);
hy1=h*y1;
fprintf('\n%0.1f %0.4f ',x,y);
y=y+hy1;
x=x+h;
end
fprintf('\n');
fprintf('\n');
disp('concentracion de CB');
x=0;xf=2;y=1;
h=0.1;
n=(xf-x)/h;
f=inline('1*x-0.5*y','x','y');
disp(' x y');
for i=1:n+1
y1=feval(f,x,y);
hy1=h*y1;
fprintf('\n%0.1f %0.4f ',x,y);
y=y+hy1;
x=x+h;
end
EJERCICIO 5.7:
Resuelva los EDO del ejercicio 5.7 utilizando el programa de Euler VBA
como dan en el ejemplo 5.7. Use los pasos de tiempo de 0.1, 0.01 y 0.001
minutos, respectivamente y compare los resultados para encontrar una
solución aceptable.
Variar los valores de h
Solución:
Variar los valores de h en el algoritmo.
fprintf('\n\n\n\t\t\t\t MÉTODO DE EULER ')
clear all; clc
syms x
syms y
x=0;xf=1;y=1;
h=0.1;
n=(xf-x)/h;
f=inline('-1*x+0.5*y','x','y');
disp('CONCENTRACION CA')
disp(' x y');
for i=1:n+1
y1=feval(f,x,y);
hy1=h*y1;
fprintf('\n%0.1f %0.4f ',x,y);
y=y+hy1;
x=x+h;
end
fprintf('\n');
fprintf('\n');
disp('concentracion de CB');
x=0;xf=2;y=1;
h=0.1;
n=(xf-x)/h;
f=inline('1*x-0.5*y','x','y');
disp(' x y');
for i=1:n+1
y1=feval(f,x,y);
hy1=h*y1;
fprintf('\n%0.1f %0.4f ',x,y);
y=y+hy1;
x=x+h;
end
EJERCICIO 5.8:
implementar el método de Runge-Kutta de segundo orden para cualquier cantidad de
ODE-IVP simultáneos en VBA. Use el ejemplo 5.7 como modelo para diseñar, codificar y
probar el programa. Resuelva el problema del ejercicio 5.7 para probar el programa.
Solución
Runge – Kutta en Matlab (Ode23)
clear all,clc
% k1=1;
% k2=0.5;
Eq=inline('[-1*x(1)+0.5*x(2);1*x(1)-0.5*x(2)]','t','x'); %ecuaciones de
diseño de reactor (se pueden cambiar para cualquier problema)
[t,x]=ode23(Eq,[0:0.1:2],[1 0]); %dt=0.1 (se puede cambiar para cualquier
problema)
disp([t,x])
plot(t,x(:,1),'-r',t,x(:,2),'-b','LineWidth',2),grid on
title('Solución de un sistema de ecuaciones reactivas','Fontsize',15);
text(0.2,0.8,'CA');text(0.8,0.1,'CB');
xlabel('Tiempo');ylabel('Concentración');
EJERCICIO 5.9:
Fermentación con Penicilina
Un modelo para un reactor discontinuo en el que la penicilina se produce por
fermentación se ha derivado de la siguiente manera (constantantinides et al.1970) para la
producción celular y la síntesis de penicilina, respectivamente:
𝑑𝑦1 𝑏1
= 𝑏1 𝑦1 − 𝑦1 2 𝑦1 (0) = 0.03
𝑑𝑡 𝑏2
𝑑𝑦2
= 𝑏3 𝑦1 𝑦2 (0) = 0.0
𝑑𝑡
Donde
𝑦1 =Concentración adimensional de la masa celular
𝑦2 =Concentración adimensional de la penicilina
𝑡 = tiempo adimensional, 0 ≤ 𝑡 ≤ 1
Los experimentos han determinado que:
𝑏1 = 13.1
𝑏2 = 0.94
𝑏3 = 1.71
a. Resolver este ODE-IVP utilizando el método de Euler y Excel (no VBA).
b. Resolver este ODE-IVP utilizando el método RK de segundo orden y Excel (no VBA)
c. Resolver este ODE-IVP usando el programa Euler VBA como se da en el ejemplo
5.7.
d. Resolver este ODE-IVP utilizando el segundo orden RK VBA programa de ejercicio
5.8.
En todos los casos, experimente con el paso del tiempo para asegurar una solución exacta.
También, gráfica 𝑦1 y 𝑦2 en función del tiempo con anotaciones apropiadas.
SOLUCION
B.-Usando el método de Runge-Kutta 2° orden
Usando el software Polymath para resolver sistemas de ecuaciones diferenciales
se obtienen los siguientes resultados:
𝑦1 = 0.939942
𝑦2 = 1.184737
Grafica
Código del programa:
o Usando el software matlab se obtienen los siguientes resultados
Código matlab:
clc, clear all
fprintf('\t********************************************\n')
fprintf('\t*\tCURSO: METODOS NUMERICOS \t*\n')
fprintf('\t*\tPROFESOR: Dr. Guillermo Evangelista B.\t*\n')
fprintf('\t*\tRESOLUCION DEL EJERCICIO 5.9 \t*\n')
fprintf('\t********************************************\n\n')
fprintf('\t Parte (a)\n')
format short
f=inline('[13.1*y(1)-
((13.1/0.94)*y(1)^2);1.71*y(1)]','t','y');
[t,y]=ode23(f,[0:0.0025:1],[0.03 0]);
tab=[t,y];
fprintf('\t\tMETODO DE RUNGE-KUTTA 2° ORDEN\n')
fprintf('\t-----------------------------------\n')
fprintf('\t\t t y1 y2\n')
fprintf('\t-----------------------------------\n')
fprintf('\t%8.4f %10.5f %10.5f\n',tab')
fprintf('\t-----------------------------------\n\n')
plot(t,y(:,1),'-r',t,y(:,2),'-b','LineWidth',1)
title('Fermentacion con penicilina')
xlabel('Tiempo, t')
ylabel('Concentración, y')
text(0.35,0.8,'y1','color','m');text(0.3,0.2,'y2','color','m'
);
legend('y1: Concentracion de la masa celular RK-2','y2:
Concentracion de la penicilina RK-2')
grid on
Usando este software se obtienen los siguientes resultados
𝑦1 = 0.93999
𝑦2 = 1.18453
Grafica
D.- Usando el método Runge-Kutta 4° orden
Código del programa:
g=inline('[13.1*y(1)-
((13.1/0.94)*y(1)^2);1.71*y(1)]','t','y');
[t,y]=ode45(g,[0:0.0025:1],[0.03 0]);
tab1=[t,y];
fprintf('\t Parte (b)\n')
fprintf('\t\tMETODO DE RUNGE-KUTTA 4° ORDEN\n')
fprintf('\t-----------------------------------\n')
fprintf('\t\t t y1 y2\n')
fprintf('\t-----------------------------------\n')
fprintf('\t%8.4f %10.5f %10.5f\n',tab1')
fprintf('\t-----------------------------------\n\n')
hold on
plot(t,y(:,1),'-g',t,y(:,2),'-m','LineWidth',1)
title('Fermentacion con penicilina')
xlabel('Tiempo, t')
ylabel('Concentración, y')
text(0.35,0.8,'y1','color','m');text(0.3,0.2,'y2','color','m'
);
legend('y1: Concentracion de la masa celular RK-2','y2:
Concentracion de la penicilina RK-2','y1: Concentracion de la
masa celular Rk-4','y2: Concentracion de la penicilina RK-4')
Usando este software se obtienen los siguientes resultados
𝑦1 = 0.93994
𝑦2 = 1.18485
Grafica
EJERCICIO 5.10:
Cinética de Bioreacciones
El modelo "Monod" para la cinética de bioreacción puede expresarse como:
𝑑𝑠 𝑘𝑠𝑥
=−
𝑑𝑡 𝑘𝑠 + 𝑠
𝑠(0) = 𝑠0 𝑥(0) = 𝑥0
𝑑𝑥 𝑘𝑠𝑥
=𝑦 − 𝑏𝑥
𝑑𝑡 𝑘𝑠 + 𝑠
Donde:
𝑠 = Concentración de sustrato limitante del crecimiento (𝑀𝐿−3 )
𝑥 = 𝑐𝑜𝑛𝑐𝑒𝑛𝑡𝑟𝑎𝑐𝑖𝑜𝑛 𝑑𝑒 𝑏𝑖𝑜𝑚𝑎𝑠𝑎 (𝑀𝐿−3 )
𝑘 = 𝑉𝑒𝑙𝑜𝑐𝑖𝑑𝑑 𝑑𝑒 𝑐𝑜𝑚𝑠𝑢𝑚𝑜 𝑚𝑎𝑥𝑖𝑚𝑎 𝑑𝑒𝑙 𝑠𝑢𝑠𝑡𝑟𝑎𝑡𝑜(𝑇 −1 ) = 5
𝑘𝑠 = 𝐶𝑜𝑛𝑠𝑡𝑎𝑛𝑡𝑒 𝑑𝑒 𝑚𝑒𝑑𝑖𝑎 𝑠𝑎𝑡𝑢𝑟𝑎𝑐𝑖ó𝑛 𝑝𝑎𝑟𝑎 𝑒𝑙 𝑐𝑟𝑒𝑐𝑖𝑚𝑖𝑒𝑛𝑡𝑜(𝑀𝐿−3 ) = 20
𝑦 = 𝑐𝑜𝑒𝑓𝑖𝑐𝑖𝑒𝑛𝑡𝑒 𝑑𝑒 𝑟𝑒𝑛𝑑𝑖𝑚𝑖𝑒𝑛𝑡𝑜(𝑀𝑀−1 ) = 0.05
𝑏 = 𝑐𝑜𝑒𝑓𝑖𝑐𝑖𝑒𝑛𝑡𝑒 𝑑𝑒 𝑑𝑒𝑐𝑎𝑖𝑚𝑖𝑒𝑛𝑡𝑜 (𝑇 −1 ) = 0.01
Las condiciones iniciales son 𝑠0 = 1000 𝑦 𝑥0 = 100
a. Resolver este ODE-IVP utilizando el método de Euler y Excel (no VBA).
b. Resolver este ODE-IVP utilizando el método RK de segundo orden y Excel (no VBA)
c. Resolver este ODE-IVP usando el programa Euler VBA como se da en el ejemplo
5.7.
d. Resolver este ODE-IVP utilizando el segundo orden RK VBA programa de ejercicio
5.8.
Un tiempo máximo apropiado para integrar estas ecuaciones debe ser determinado por
experimentación. También, experimente con el paso del tiempo para garantizar buenos
resultados. Gráfique 𝑠 y 𝑥 en función del tiempo con las anotaciones apropiadas (hacer
esto con datos para sólo un paso de tiempo apropiado).
SOLUCION
Usando el software Polymath para resolver sistemas de ecuaciones diferenciales
se obtienen los siguientes resultados:
Grafica
Para (a)
Usando el software MATLAB para resolver sistemas de ecuaciones diferenciales se
obtienen los siguientes resultados:
Grafica
Código del programa:
function [t,y,x]=euler(f1,f2,t0,tf,y0,x0,n)
h = (tf-t0)/n;
t = t0:h:tf;
x(1)=x0;
y(1)=y0;
for i=1:n
x(i+1)=x(i)+h*f1(x(i),y(i));
y(i+1)=y(i)+h*f2(x(i),y(i));
end
f1=@(x,y) 0.05*x*y/(20+y);
f2=@(x,y) -5*y*x/(20+y);
[t,y,x]=euler(f1,f2,0,2,1000,100,15);
hold on
plot(t,x,'b')
plot(t,y,'r')
xlabel('t')
ylabel('x,s')
axis([0 2 0 1000])
title('sistema de ecuaciones diferenciales')
hold off
grid
tab=[t',y',x'];
fprintf('\t\tMETODO DE EULER\n')
fprintf('\t-----------------------------------\n')
fprintf('\t\t t s x\n')
fprintf('\t-----------------------------------\n')
fprintf('\t%8.4f %10.5f %10.5f\n',tab')
fprintf('\nse observa en la grafica que para un tiempo aproximadamente de
t=2h nos garantiza\n')
fprintf('una buena conversion de la biorreaccion.\n')
Para (b)
Usando el software MATLAB para resolver sistemas de ecuaciones diferenciales se
obtienen los siguientes resultados:
Grafica
Código del programa
y0(1)=100;
y0(2)=1000;
tspan=[0 2];
fg=@(t,y) [0.05*y(1)*y(2)/(20+y(2));-5*y(2)*y(1)/(20+y(2))];
[t,y]=ode23(fg,tspan,y0);
tab=[t,y];
fprintf('\t\tMETODO DE RUNGE-KUTTA 2° ORDEN\n')
fprintf('\t-----------------------------------\n')
fprintf('\t\t t x s\n')
fprintf('\t-----------------------------------\n')
fprintf('\t%8.4f %10.5f %10.5f\n',tab')
fprintf('\t-----------------------------------\n\n')
plot(t,y(:,1),'-r',t,y(:,2),'-b','LineWidth',1)
title('bioreaccion')
xlabel('Tiempo, t')
ylabel('Concentración,x, s')
fprintf('\nse observa en la grafica que para un tiempo aproximadamente de
t=2h nos garantiza\n')
fprintf('una buena conversion de la biorreaccion.\n')
EJERCICIO 5.11:
implementar el método de Runge-Kutta de cuarto orden para cualquier cantidad de
ODE-IVP simultáneos en VBA.
Lo único que un usuario debe cambiar para usar el programa es Subrountine para
calcular las funciones del "lado derecho" para las EDO. El usuario también debe
especificar los datos de entrada requeridos. Una interfaz de usuario sugerida asociada
con la hoja de cálculo de Excel que interactúa con el programa VBA es como se muestra
en el ejemplo 5.7. (The only thing a user must change to use the program is the
Subrountine to calculate the "right-hand side" functions for the EDOs. The user must
also specify the required input data. A suggested user interface associated with the
Excel spreadsheet that interacts with the VBA program is as shown in example 5.7.)
A continuación, se dan algunos consejos para ayudar a preparar el programa VBA.
Recuerde el algoritmo para usar el método de cuarto orden [Link] con solo un
EDO:
k1 = hf( yn, tn)
𝒌𝟐 𝒉
k2 = hf( yn + , t n +𝟐 )
𝟐
𝒌𝟐 𝒉
k3 = hf( yn + , t n +𝟐 )
𝟐
k4 = hf( yn + k3, tn + h )
𝟏
yn+1 = yn + 𝟔 [𝒌𝟏 + 𝟐𝒌𝟐 + 𝟐𝒌𝟑 + 𝒌𝟒 ]
Para generalizar esto a N EDO simultáneas, haga que cada k, y, y f una matriz con
subíndices de 1 a N. Aquí están las declaraciones ReDim sugeridas:
ReDim k1 (N) como doble
ReDim k2 (N) como doble
ReDim k3 (N) como doble
ReDim k4 (N) como doble
ReDim y (N) como doble
ReDim f (N) como doble
ReDim z (N) como doble
Use una subrutina (como Sub FCalc del Ejemplo 5.7) que tenga como entrada
la hora actual, la matriz y, y el número de ecuaciones (N), y regresa
N funciones del lado derecho (f). Esta subrutina se debe llamar cuatro veces
por paso de tiempo Tenga en cuenta que y, k1, k2, k3 y k4 son vectores de N-longitud.
1. Una vez en el punto base f( yn, tn)
𝒌 𝒉
2. Una primera vez en el punto medio f( yn + 𝟐𝟏 , tn +𝟐 )
𝒌 𝒉
3. Una segunda vez en el punto medio f( yn + 𝟐𝟐 , tn +𝟐 )
4. Y finalmente en el punto final f( yn + k3, tn + h )
Se sugiere definir otra matriz (llámala z) y usar esto para el primer
argumento en las llamadas de subrutina a mitad de camino y puntos finales. Además,
defina una segunda variable de tiempo (llámalo Ttemp). Por ejemplo, para la primera
llamada al punto medio, usa la secuencia.
z(i) = y(i) + k1(i)/2 i = 1, 2,..., N
Ttemp = Time + h/2
Call FCalc(Ttemp, z(), N, f())
k2(i) = h*f(i) i = 1, 2,..., N
Se recomienda encarecidamente esbozar la lógica del código y "ejecutar
a través de "esta lógica cuidadosamente antes de comenzar realmente a escribir el
código.
Finalmente, pruebe el programa VBA resultante resolviendo el ODE-IVP de
Ejercicio 5.10. Experimente con el paso del tiempo para asegurarse de que un aceptable
la solución está siendo generada. Tenga cuidado con los resultados para un paso de
tiempo de 0.1.
Algoritmo de método Range – Kutta de 4to orden
xf=input(‘ ‘)
h= input(‘ ‘)
xo= input(‘ ‘);
yo= input(‘ ‘); %condicion inicial
i= input(‘ ‘);
f=inline('[ingresar
ecuación]','x','y');
x1=zeros(1,length(N));
x1(1,1)=xo;
y1=zeros(1,length(N));
y1(1,1)=yo;
fprintf('\n\n\t\t\t\t\tpara dt = 0.05\n')
fprintf('\t\t\t\t------------------------------\n')
fprintf('\t\t\t \ti \t t(i) \tH(i) \n')
fprintf('\t\t\t\t------------------------------\n')
fprintf('\t\t\t \t%2d \t%6.4f \t%8.6f\n',i,xo,yo)
N=(xf-xo)/h;
while i<N
k1=f(xo,yo);
k2=f(xo+h/2,yo+h*k1/2);
k3=f(xo+h/2,yo+h*k2/2);
k4=f(xo+h,yo+h*k3);
yo=yo+h/6*(k1+2*k2+2*k3+k4);
xo=xo+h;
i=i+1;
fprintf('\t\t\t \t%2d \t%6.4f \t%8.6f\n',i,xo,yo)
x1(1,i+1)=xo;
y1(1,i+1)=yo;
end
Algoritmo Usando Runge – Kutta de 4to Grado en MatLab (Ode45)
y0(1)=100; %valores iniciales
y0(2)=1000;
tspan=[0 2];
fg=@(t,y) [0.05*y(1)*y(2)/(20+y(2));-5*y(2)*y(1)/(20+y(2))];%ingreso de
ecuaciones
[t,y]=ode45(fg,tspan,y0); %Runge-Kutta de 4to orden (ode45)
tab=[t,y];
fprintf('\t\tMETODO DE RUNGE-KUTTA 4° ORDEN\n')
fprintf('\t-----------------------------------\n')
fprintf('\t\t t x s\n')
fprintf('\t-----------------------------------\n')
fprintf('\t%8.4f %10.5f %10.5f\n',tab')
fprintf('\t-----------------------------------\n\n')
plot(t,y(:,1),'-r',t,y(:,2),'-b','LineWidth',1)
title('bioreaccion')
xlabel('Tiempo, t')
ylabel('Concentración,x, s')
fprintf('\nse observa en la grafica que para un tiempo aproximadamente de
t=2h nos garantiza\n')
fprintf('una buena conversion de la biorreaccion.\n')
BIBLIOGRAFIA
Chapra, S. y Canale, R. 2015. Métodos numéricos para ingenieros. Séptima
Edició[Link]-Hill/Interamericana Editores, S.A. de C.V. México
Nieves, A. y Domínguez, F. 2012. Métodos Numéricos Aplicados a la Ingeniería.
Cuarta Edición. Grupo Editorial Patria, S.A. de C. V. México.
