SOLUCION_MODELAMIENTO

Realizado por: JUVENAL TORDOCILLO PUCHUC Lic. en Física

Modelamiento Procesos Ambientales

Mg. en Física con mención en geofísica

INTRODUCCION El curso de modelamiento de procesos ambientales, se ha utilizado el software Matlab para simular distintos procesos ambientales tales como contaminación ríos, lago y atmosfera. Como parte del desarrollo del curso y del nivel alcanzado en el manejo del Matlab para procesos ambientales se ha desarrollado 5 ejercicios con respectivos script adjunto a este archivo donde en cada solución se ha puesto énfasis que cada transformación de datos sea lo más automático posible y para el desarrollo de sistema de ecuaciones ordinarias se procedió mediante el uso del método de Euler, si bien es un método básico pero permite encontrar resultados satisfactorios para incrementos pequeños y siempre en cuando la precisión de encontrar los datos no sean tan relevantes, como es el caso este trabajo.

-1-




PROBLEMA PROPUESTO: DE AJUSTE DE DATOS EXPERIMENTALES

Enunciado 1: En la ingeniería de abastecimiento de aguas, el tamaño del reservorio depende de la estimación exacta del flujo de agua en el río del cual se toma. Para algunos rios es difícil obtener registros históricos de muchos años atrás de tales datos de flujo. Por el contrario, datos meteorológicos sobre precipitación de muchos años atrás están a menudo disponibles. Por tanto, con frecuencia es útil determinar una relación entre el flujo y precipitación. Esta relación se puede entonces usar para estimar flujos por años pero solo cuando se hicieron dichas mediciones de precipitación. Para un río que se va a encauzar a un dique, se tienen los siguientes datos:

Precipitación, cm Flujo, m3/s 88.9

114.7

101.6

172.0

104.1

152.9

139.7

269.0

132.1

206.4

94.0

161.4

116.8

175.8

121.9

239.0

99.1

130.0

a) Grafique los datos b) Ajuste a una línea recta. Sobreponga la línea de su gráfica c) Use la mejor línea de ajuste para predecir el flujo de agua anual si la precipitación es de 120 cm.

SOLUCION a) Grafique los datos b) Ajuste a una línea recta. Sobreponga la línea de su gráfica

El procedimiento que permite leer los datos externos ejercicio_1 txt es el comando load.

-2-




PASO 1: LECTURA DE DATOS clc clear all % procedimiento de lectura y graficado de los datos % este comando permite ordenar de forma ascendente A=load('ejercicio_2.txt');

PASO 2: Ordenar en orden ascendente se utilizó el siguiente procedimiento [B,K]=sort(A(:,2)); B = [A(K) B]; x=B(:,1); y=B(:,2); plot(x,y,'bs','linewidth',1) grid on xlabel('PRECIPITACION (cm)') ylabel('FLUJO (m^3/s)') xmin=min(x); xmax=max(x);

PASO 3: Ajuste de datos con polinomio de primer orden y mostrar gráfico % procedimiento de ajuste de los datos %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% figure () % gafico mostrado para ajuste de orden 1 n=input('Ingrese el grado de polinomio 1:'); p1=polyfit(x,y,n); x1=linspace(xmin,xmax,50); y1=polyval(p1,x1); plot(x,y,'bp',x1,y1,'ro-','linewidth',1) legend('datos experimentales','ajuste orden 1') xlabel('PRECIPITACION (cm)') ylabel('FLUJO (m^3/s)') grid on

PASO 4: Ajuste de datos con polinomio de segundo orden y mostrar gráfico -3-




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% figure () % gafico mostrado para ajuste de orden 2 n=input('Ingrese el grado de polinomio 2:'); p2=polyfit(x,y,n); x2=linspace(xmin,xmax,50); y2=polyval(p2,x2); plot(x,y,'bp',x2,y2,'r*-') legend('datos experimentales','ajuste orden 2') grid minor xlabel('PRECIPITACION (cm)') ylabel('FLUJO (m^3/s)')

PASO 5: Ajuste de datos con polinomio de cuarto orden y mostrar gráfico figure () % gafico mostrado para ajuste de orden 4 n=input('Ingrese el grado de polinomio 4:'); p3=polyfit(x,y,n); x3=linspace(xmin,xmax,50); y3=polyval(p3,x3); plot(x,y,'bp',x3,y3,'r*-','linewidth',1) legend('datos experimentales','ajuste orden 4') grid minor xlabel('PRECIPITACION (cm)') ylabel('FLUJO (m^3/s)') %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

PASO 6: Muestra de gráfico figure () % gafico mostrado todos los ajustes plot(x,y,'bp',x1,y1,'r*-',x2,y2,'b+-',x3,y3,'gp-') legend('datos experimentales','ajuste orden 1','ajuste orden 2','ajuste orden 4') grid minor xlabel('PRECIPITACION (cm)') ylabel('FLUJO (m^3/s)')

-4-




c) Use la mejor línea de ajuste para predecir el flujo de agua anual si la precipitación es de 120 cm.

PASO 7: Calculo de datos para el problema %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% p11=polyval(p1,124); fprintf('AJUSTE CON POLINOMIO DE GRADO 1:') fprintf('\t\t %10.2f\n\n',p11) p22=polyval(p2,124); fprintf('AJUSTE CON POLINOMIO DE GRADO 2:') fprintf('\t\t %10.2f\n\n',p22) p33=polyval(p3,124); fprintf('AJUSTE CON POLINOMIO DE GRADO 4:') fprintf('\t\t %10.2f\n\n',p33)

-5-




Enunciado 2: Tres organismos portadores de enfermedades decaen de manera exponencial en las aguas de un lago de acuerdo con el siguiente modelo: p(t) = Ae-1.5t + Be- 0.3t + Ce-0.05t Calcule la población inicial de cada organismo (A, B y C) dadas las siguientes mediciones: t, HR

p(T)

0.5

7

1

5.2

2

3.8

3

3.2

4

2.5

5

2.1

6

1.8

7

1.5

8

1.2

9

1.1

SOLUCION El procedimiento que permite leer los datos externos ejercicio_2 txt es el comando load. PASO 1: LECTURA DE DATOS clc clear all clf % procedimiento de lectura y graficado de los datos % este comando permite ordenar de forma ascendente A=load('ejercicio_2.txt');

PASO 2: MOSTRA GRAFICO x=A(:,1); y=A(:,2);

-6-




plot(x,y,'ro') grid on xlabel('TIEMPO(HORAS)') ylabel('POBLACION (MICRO-ORGAMISMOS)')

PASO 3: Ajuste de datos exponenciales y mostrar gráfico %pause % procedimiento de ajuste de los datos hold on f2 = fit(x,y,'exp2') plot(f2,x,y) xlabel('TIEMPO(HORAS)') ylabel('POBLACION (MICRO-ORGAMISMOS)') %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

PASO 4: Transformación de datos a lineal %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% figure () % gafico mostrado para ajuste de orden 3 y1=log(y); p1=polyfit(x,y1,1);

-7-




xmin=min(x); xmax=max(x); xx=linspace(xmin,xmax,50); yy=polyval(p1,xx); plot(x,y1,'ro','linewidth',2) hold on plot(xx,yy,'k--') legend('datos experimentales','ajuste orden 1') xlabel('TIEMPO(HORAS)') ylabel('LN(POBLACION (MICRO-ORGAMISMOS))') grid on

PASO 5: Encuentre el coeficiente correlación %coeficente de correlacion lineal CC=corrcoef(x,y1); rxy=CC(1,2); rr2=rxy*rxy; fprintf('COEFICIENTE DE CORRELACION LINEAL :') fprintf('\t\t %10.2f\n\n',rr2)

PASO 6: Encuentre la raíz y de corte en las abcias y la ordenada %--------------------------------------------%Encuentra la raiz en la abcisa y11=0; f = @(xx) polyval(p1,xx) - y11; raiz=fzero(f,5); fprintf('LA INTERSECCION EN EL EJE ABCISAS (TIEMPO):') fprintf('\t\t %10.2f\n\n',raiz) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% p11=polyval(p1,0); fprintf('LA INTERSECCION CON EL EJE DE LA ORDENADA:')

-8-




fprintf('\t\t %10.2f\n\n',p11)

PASO 7: Calculo de A, B, C A1='exp(1.85)=A+B+C'; B1='1=A*exp(-1.5*8.9)+B*exp(-0.3*8.9)+A*exp(-0.05*8.9)' C1='7=A*exp(-1.5*0.5)+B*exp(-0.3*0.5)+A*exp(-0.05*0.5)' [A,B,C]=solve(A1,B1,C1,'A,B,C'); disp('MICROORGANISMO A:'),A disp('MICROORGANISMO B:'),B disp('MICROORGANISMO C:'),C

-9-




Enunciado 3: Modelos Predador-Presa con Migración En este caso, asumimos que las presas constituyen una especie migratoria. Para ello se añade al modelo LV un término que refleja la emigración e inmigración periódica. Dicho término puede ser una función sinusoidal con amplitud M y periodo w:

dx  Ax  Bxy  Msen(wt ) dt dy  Cy  Dxy dt

a. Sea que las presas migran todos los años entre este santuario y otra zona en la que crían, retornando posteriormente. Analizar el efecto de esta migración sobre el comportamiento del modelo. SOLUCION PASO 1: LECTURA DE DATOS clc clear all % Entrada de datos y declarcion de parametros h=0.01; ti=input('Ingrese el tiempo inicial: '); tf=input('Ingrese el tiempo final: '); N=(tf-ti)/h; A=1.0; B=0.5; C=0.75; D=0.25; X(1)=1.5; Y(1)=2; W=0.5; M=0.01; t=ti:h:tf;

PASO 2: TECNICA NUMERICA DE EULER % sE UTILIZARA EL METODO DE EULER PARA APROXIMACION for i=1:N dx=A*X(i)-B*X(i)*Y(i)-M*sin(W*t(i)); dy=-C*Y(i)+D*X(i)*Y(i); X(i+1)=X(i)+h*dx; Y(i+1)=Y(i)+h*dy; end

- 10 -




PASO 4: Muestra del gráfico y su interpretación plot(t,X,'r-') ylabel('PRESA') xlabel('tiempo') title ('MODELO LOGISTICO LV') grid on

PASO 4: Muestra del gráfico y su interpretación figure() plot(t,Y,'k-') ylabel('PREDADOR') xlabel('tiempo') title ('MODELO LOGISTICO LV') grid on

- 11 -




PASO 5: Muestra del gráfico y su interpretación figure () plot(t,X,'r-', t,Y,'k-') ylabel('PRESA y PRESA') xlabel('tiempo') title ('MODELO LOGISTICO LV') legend('PRESA','PREDADOR') grid on

PASO 6: Muestra del gráfico y su interpretación figure () plot(t,X,'r-', t,Y,'k-') ylabel('PRESA y PRESA') xlabel('tiempo') title ('MODELO LOGISTICO LV') legend('PRESA','PREDADOR') grid on

- 12 -




Enunciado 4 : TRES CSTRs EN SERIE Las ecuaciones que describen a tres reactores CSTR isotérmicos en serie son dC A1 1 = (CA0 – CA1) – k CA1  dt dC A2 1 = (CA1 – CA2) – k CA2  dt dC A3 1 = (CA2 – CA3) – k CA3  dt Las condiciones iniciales son CA1(0) = 0,4 kg.mol de componente A/m3, CA2(0) = 0,2 kg.mol de componente A/m3, y CA3(0) = 0,1 kg.mol de componente A/m3. La función impulsora es CA0 . Asumiremos que en el tiempo cero CA0 se ajusta a 1,8 kg.mol de A/m3 y se mantiene constante. El parámetro  = V/F (tiempo de residencia), es ajustado igual a 2 min y el valor de k es 0,5 min–1. La perturbación es el cambio en escalón de la concentración de la alimentación en el tiempo igual a cero desde 0.8 hasta 1.8 kg mol de A/m3. El tiempo está en minutos Similar al anterior y basado en el método de EULER PASO 1: Declaración de parámetros clc clear all h=0.01; ti=input('Ingrese el tiempo inicial: '); tf=input('Ingrese el tiempo final: '); N=(tf-ti)/h; Ca1(1)=0.4; Ca2(1)=0.2; Ca3(1)=0.1; Cao=1.8; k=0.5; tau=2; t=ti:h:tf;

PASO 2: Método de EULER % sE UTILIZARA EL METODO DE EULER PARA APROXIMACION for i=1:N dCa1=(Cao-Ca1(i))/tau -k*Ca1(i); dCa2=(Ca1(i)-Ca2(i))/tau -k*Ca2(i); dCa3=(Ca2(i)-Ca3(i))/tau -k*Ca3(i); Ca1(i+1)=Ca1(i)+h*dCa1; Ca2(i+1)=Ca2(i)+h*dCa2; Ca3(i+1)=Ca3(i)+h*dCa3; End

- 13 -




PASO 3: GRAFICOS UTILIZADO SUBPLOT figure subplot(2,2,1) plot(t,Ca1,'r') title('Subplot 1: Ca1') ylabel('Ca1 (kg.mol)') xlabel('tiempo (min.)') grid on subplot(2,2,2) plot(t,Ca2,'r') title('Subplot 2: Ca2') ylabel('Ca2 (kg.mol)') xlabel('tiempo (min.)') grid on subplot(2,2,3) plot(t,Ca3,'r') title('Subplot 3: Ca3') ylabel('Ca1 (kg.mol)') xlabel('tiempo (min.)') grid on subplot(2,2,4) plot(Ca1,Ca2,'k') title('Subplot 4: Ca1 vs Ca2') ylabel('Ca1 (kg.mol)') xlabel('Ca2 (kg.mol)') grid on

PASO 4: RESULTADOS

PASO 5: RESULTADOS figure () subplot(2,1,1); plot(Ca2,Ca3,'k') title('Subplot 4: Ca2 vs Ca3') ylabel('Ca2 (kg.mol)') xlabel('Ca3 (kg.mol)') grid on

- 14 -




subplot(2,1,2); plot(Ca3,Ca1,'k') title('Subplot 4: Ca1 vs Ca3') ylabel('Ca1 (kg.mol)') xlabel('Ca3 (kg.mol)') grid on

Enunciado 5 : PROGRAMACIÓN LINEAL “MODELAMIENTO Y SIMULACIÓN DE UN CASO DE TRATAMIENTO DE RESIDUOS CONTAMINANTES” La distribución bidimensional de la concentración de contaminantes en un canal se puede describir con

c(x, y) = 7,9 + 0,13x + 0,21 y – 0,05 x2 – 0,016 y2 – 0,007 xy Determine la localización exacta de La concentración pico dada la función y con el conocimiento de que el pico cae dentro de las fronteras –10  x  10 y 0  y  20. PASO 1: DECLARACION DE DATOS DE X e Y ADEMAS DE LA ECUACION CORRESPONDIENTE clc clear all x=linspace(-10,10,50); y=linspace(0,20,50); [x1 y1]=meshgrid(x,y); C=7.9+0.13*x1+0.21*y1-0.05*x1.^2-0.016*y1.^2-0.007*x1.*y1;

- 15 -




PASO 2: GRAFICO EN 3D surf(x1,y1,C) title('CONCENTRACION EN UN CANAL BIDIMENSIONAL') ylabel('X (m)') xlabel('Y (m)') shading interp

PASO 3: ENCUENTRA EL MAXIMO Y LAS COORDENADAS QUE CORRESPONDE xmax=max(max(C)); xmin=min(min(C)); [fila col]=find(C==max(max((C)))) fprintf('EL VALOR MAXIMO DE LA CONCENTRACIÓN ES :') fprintf('\t\t %10.2f\n\n',xmax) fprintf('en el eje y es :') fprintf('\t\t %10.2f\n\n',fila) fprintf('en el eje x es :') fprintf('\t\t %10.2f\n\n',col) hc = colorbar; set(get(hc,'title'),'string','Concentracion','FontS',11,'FontW','Bold' ) set(hc,'YTick',[xmin:xmax],'YTickLabel',[xmin:xmax],'FontS',10)

CONCLUSION 

Se ha utilizado el software matlab para implementar los modelos con enfoque al comportamiento de contaminación ambiental.



Se realizado todo el proceso en base a la implementación de script que puede ser modificable según los parámetros se requiera modificar. - 16 -

SOLUCION_MODELAMIENTO

Recommend Documents