- Matlab -
TALLER 4
1) Animación
% Movimiento Browniano
n = 20;
s = .02;
x = rand(n,1)-0.5;
y = rand(n,1)-0.5;
h = plot(x,y,'.');
axis([-1 1 -1 1])
axis square
grid off
set(h,'EraseMode','xor','MarkerSize',18)
while 1
drawnow
x = x + s*randn(n,1);
y = y + s*randn(n,1);
set(h,'XData',x,'YData',y)
end
2) Grabación de animaciones en formato avi.
clear; clc;
fig=figure;
set(fig,'DoubleBuffer','on');
figure(1);
set(gca,'xlim',[-80 80],'ylim',[-80 80],'NextPlot','replace','Visible','off')
mov = avifile('EjmAVI1.avi')
x = -pi:.1:pi;
radio = 0:length(x);
for k=1:length(x)
h = patch(sin(x)*radio(k),cos(x)*radio(k),[abs(cos(x(k))) 0 0]);
set(h,'EraseMode','xor');
F = getframe(gca);
mov = addframe(mov,F);
end
mov = close(mov);
3) Minimización no lineal en un intervalo
a) x = fminbnd(@sin,0,2*pi)
b) x = fminbnd(@FMin,0,5)
function f = FMin(x)
f = (x-3).^2 - 1;
- Matlab -
4) Minimización con restricciones
Hallar los valores de x que minimizan f ( X ) = − X 1 X 2 X 3
Sujeto a las siguientes restricciones:
− X1 − 2 X 2 − 2 X 3 ≤ 0
X 1 + 2 X 2 + 2 X 3 ≤ 72
x0 = [10; 10; 10]; % Supuestos
A = [-1 -2 -2; 1 2 2];
b = [0 72]’;
options=optimset('Diagnostics','off','Display','Iter','LineSearchType','cubicpoly','Hessian'
,'off');
[x,fval] = fmincon(@FMinR,x0,A,b,[],[],[],[],[],options)
function f = FMinR(x)
f = -x(1) * x(2) * x(3);
5) Hallar el vector x que minimice:
∑ ( 2 + 2k − e − e kx2 )
10
kx1
k =1
iniciando en el punto x = [0.3, 0.4].
x0 = [0.3 0.4] % puntos iniciales
[x,resnorm] = lsqnonlin(@FMinMC,x0) % Invocar al optimizador
function F = FMinMC (x)
k = 1:10;
F = 2 + 2*k-exp(k*x(1))-exp(k*x(2));
6) PROCESO BIOLÓGICO EN UN BIOREACTOR
CONCEPTOS DEMOSTRADOS
Simulación de reacciones químicas o biológicas en un proceso por lotes que puede tratar
con velocidades muy altas de reacción con muy bajas concentraciones de reactantes.
METODOS NUMERICOS USADOS
Solución de sistemas de ecuaciones diferenciales ordinarias de primer orden que llegan
a ser difíciles durante el curso de la integración.
ENUNCIADO DEL PROBLEMA
Un proceso biológico involucra el crecimiento desde el sustrato. Los balances de
materia en este proceso batch producen:
- Matlab -
dB BS dS BS
= k1 = −0.75k1
dt k2 + S dt k2 + S
donde B y S son las respectivas concentraciones de biomasa y sustrato. Las constantes
cinéticas de la reacción son k1=0.3 y k2=10-6 en unidades consistentes.
(a) Solucionar este conjunto de ecuaciones diferenciales empezando de t=0, donde
S=5.0 y B=0.05 hasta un tiempo final dado por tf=20. Asumir unidades consistentes.
(b) Plotear S y B con las condiciones para la parte (a).
function ReactorStiff();
clc; clf
k1=0.3; k2=1E-6; t0=0; B0=0.05; S0=5.0; tf=20;
Param=[k1 k2];
options = odeset('RelTol',1e-7,'AbsTol',1e-7);
%x(1)=B, x(2)=S
[T,Y] = ode45(@fStiff,[t0 tf],[B0 S0],options,Param);
% [T,Y] = ode15s(@fStiff,[t0 tf],[B0 S0],options,Param); %STIFF
figure(1)
plot(T,Y,'-')
title('Crecimiento de Biomasa desde el sustrato')
h = legend(['Biomasa(B)'],['Sustrato(S)'],4);
xlabel('Tiempo')
ylabel('Concentración de Biomasa (B) , Sustrato (S)')
function dx = fStiff(t,x,Param)
k1=Param(1); k2=Param(2);
dx = zeros(2,1);
%x(1)=B, x(2)=S
B=x(1); S=x(2);
dB_dt = k1*B*S/(k2+S);
dS_dt = -0.75*k1*B*S/(k2+S);
dx(1)=dB_dt; dx(2)=dS_dt;
7) DIFUSIÓN Y REACCIÓN QUÍMICA IRREVERSIBLE Y SIMULTÁNEA
Método del disparo para solucionar ecuaciones diferenciales con valores de frontera
típicamente usados en fenómenos de transporte y cinética de reacciones.
CONCEPTOS DEMOSTRADOS
Método para solucionar ecuaciones diferenciales ordinarias con valores de frontera de
dos puntos típicamente usado en fenómenos de transporte y cinética de reacciones.
- Matlab -
MÉTODOS NUMERICOS USADOS
Conversión de una ecuación diferencial de segundo orden en un sistema de dos
ecuaciones diferenciales de primer orden, un método de disparo para solucionar
ecuaciones diferenciales ordinarias (EDO), y usar el método de secante para solucionar
problemas de valor de frontera de dos puntos.
ENUNCIADO DEL PROBLEMA
La difusión y reacciones químicas irreversibles de primer orden simultáneas en una fase
simple conteniendo solo un reactante A y un producto B resulta en una ecuación
diferencial dada por
d 2CA k
2
= CA (1)
dz DAB
donde C A es la concentración del reactante A (kg-mol/mol3), z es la variable distancia
(m), k es la constante de reacción homogénea (s-1), y DAB es el coeficiente de difusión
binaria (m2/s). Una geometría típica para la ecuación anterior es de una capa
unidimensional que tiene su superficie expuesta a una concentración y no permite la
difusión a través de su superficie inferior. Así las condiciones inicial y de frontera son:
C A = C A0 para z = 0 (3)
dC A
= 0 para z = L (4)
dz
donde para C A 0 es la concentración constante en la superficie ( z = 0 ) y no hay
transporte a través de la superficie inferior (para z = L ) es decir la derivada es cero.
La ecuación diferencial tiene una solución analítica dada por
cosh ⎡⎣ L k / DAB (1 − z / L ) ⎤⎦
C A = C A0 (2)
(
cosh L k / DAB )
a) Solucionar numéricamente la ecuación (1) con las condiciones de frontera (2) y (4)
para C A 0 =0.2 kg-mol/m3, k =10-3 s-1, DAB =1.2x10-9, y L=10-3 m. Esta solución podría
usar la técnica del disparo y emplear el método de Newton-Raphson para la
convergencia a la condición de frontera dada por la ecuación (3).
b) Comparar los perfiles de concentración sobre la película como lo predice la solución
analítica de la ecuación (4).
SOLUCIÓN
function Reaction();
clc; clf
k=0.001; DAB=1.2E-9; L=1E-3; CA0=0.2;
LkDAB=L*(k/DAB)^0.5;
- Matlab -
Param=[k DAB 1];
% dCA_dz=0;%condicion en z=L debe ser 0
% Metodo de NR para satisfacer condicion de dy/dz=0 en z=L; y=dCA_dz
dh=1E-6;
dCA_dz0i=10*(-140);
dCA_dz0f=-160;
options = odeset('RelTol',1e-6,'AbsTol',1e-6);
while abs(dCA_dz0f-dCA_dz0i)>=1E-8
dCA_dz0i=dCA_dz0f;
[T,Y] = ode45(@fReaction,[0 L],[CA0 dCA_dz0i],options,Param);
fdCA_dzLCalc=Y(size(Y,1),2);
[T,Y] = ode45(@fReaction,[0 L],[CA0 (dCA_dz0i+dh)],options,Param);
f1dCA_dzLCalc=(Y(size(Y,1),2)-fdCA_dzLCalc)/dh;
dCA_dz0f=dCA_dz0i-fdCA_dzLCalc/f1dCA_dzLCalc;
end
dCA_dz0=Y(1,2);
[T,Y] = ode45(@fReaction,[0 L],[CA0 dCA_dz0],options,Param);
dCA_dzLCalc=Y(size(Y,1),2); %debe ser 0
figure(1)
subplot(2,1,1)
plot(T,Y(:,1),'r-')
title('Reacción')
xlabel('Longitud (m)')
ylabel('C_A (kgmol/mol^3), orden 1')
hold on
LZ=[0:1E-5:L];
CA=CA0*cosh(LkDAB*(1-LZ./L))./cosh(LkDAB);
plot(LZ,CA,':')
legend(['C_A (EDO)'],['C_A (analítico)'],1);
subplot(2,1,2)
plot(T,Y(:,2),'r-')
xlabel('Longitud (m)')
ylabel('dC_A/dz (kgmol/mol^3.m)')
hold on
LZ=[0:1E-5:L];
dCA_dZ=(-1/L)*LkDAB*CA0*sinh(LkDAB*(1-LZ./L))./cosh(LkDAB);
plot(LZ,dCA_dZ,':')
legend(['dC_A/dz (EDO)'],['dC_A/dz (analítico)'],2);
hold off
function dx = fReaction(t,x,Param)
k=Param(1); DAB=Param(2); Caso=Param(3);
- Matlab -
dx = zeros(2,1);
CA=x(1); y=x(2);
dCA_dz = y;
dy_z = k/DAB*CA;
dx(1)=dCA_dz;
dx(2)=dy_z;
8) PROBLEMAS PROPUESTOS
8.1 Hallar las raíces de la función f ( X ) = X 2 − 0.75sin ( X ) cos ( X ) en un rango
particular de X usando por lo menos dos métodos.
8.2 Un fluido de densidad 4 kg/L ingresa a un tanque de 120 L a un flujo de 8 L/s. El
flujo de salida es 3 L/s y el volumen inicial de 50 L. Determinar la variación de
volumen en los 100 primeros segundos usando la función ode45.
8.3 Resolver el siguiente sistema de ecuaciones en los 60 primeros segundos:
dX
= σ (−X + Y )
dt
dY
= ( λ − Z ) X − Y donde σ =10; λ =24; b =2;
dt
dZ
= XY − bZ
dt
En t=0: X=Y=Z=0.1
Plotear X vs Y vs Z en figure(1), y X vs Y, X vs Z, Y vs Z y X,Y,Z vs el tiempo en
figure(2).
Variar en forma ínfima a los 3 parámetros del sistema de ecuaciones y discutir los
resultados.
d3y ⎛ d 2 y 1⎞ ⎛ dy ⎞
− 2 x sin ⎜ ⎟ + exp ⎜ ⎟ − 3 y sin ( y ) = −t cos ( t )
3
8.4 Resolver la ecuación 4
dt ⎝ dt t ⎠ ⎝ dt ⎠
Las condiciones son y(0)=1, y’(0)=-π/2, y’’(0)=0
Usar la función ode45.
9) Leer e Imprimir Datos
9.1 Uso de INPUT: R=input('Ingrese el radio:’)
9.2 Uso de FPRINTF: fprintf('Un circulo unitario tiene una circunferencia
de %12.8f\n',2*pi)
