• Author / Uploaded
• Luis Castro

Citation preview

Solución Numérica de Ecuaciones Diferenciales Método de Euler Preparado por M.C. Luis E. Castro-Solís

¿Porqué estudiar el método de Euler? Por su sencillez, el método de Euler resulta muy adecuado para iniciar el estudio de los métodos numéricos para solución de EDO's, siendo su programación en computadora bastante abordable por estudiantes de ingeniería de nivel pre-grado. Los métodos de Euler tienen tres modalidades a) Euler hacia adelante, b) Euler hacia atras y c) Euler modificado, dependiendo si se derivan de la aproximación por diferencias hacia adelante, hacia atras o centrales respectivamente. En este documento se estudiara el caso (a) Euler hacia adelante. Problema a resolver Se desea resolver el problema de valor inicial dy ---- = F (x, y) sujeto a y(x0 ) = y0 dx en donde F(x, y) es una función que depende de "x" y de "y" (la ecuación diferencial) y "y (x0 ) = yo" es la condición incial; la solución del problema es la función y(x), lo cual significa que la condición inicial expresa lo siguiente: Cuando la variable independiente x = x0 , la variable dependiente y = y0. En resumen: Conocemos: i) La ecuación diferencial; o sea la derivada de la función desconocida y(x) en cualquier punto equivalente a F(x, y). ii) El punto inicial por el que pasa la solución y(x); o sea (x0 , y0 ). gráficamente:

Método de Euler hacia adelante

El método de Euler consiste en calcular en forma tabular aproximaciones a los puntos por los que pasa la solución y(x), usando la información del problema. Los puntos de la solución aproximada están equiespaciados (a la misma distancia horizontal unos de otros) a cada H unidades. H recibe el nombre de tamaño de paso. El procedimiento para calcular la altura aproximada y1 de la solución en x1 es: De la interpretación geométrica de la primera derivada, la pendiente de la recta tangente (se muestra en rojo en la gráfica de abajo) a la curva solución desconocida y(x) en el punto (x0 , y0 ) viene siendo la propia ecuación diferencial F(x,y) evaluada en (x0 , y0 ). Esto es, por cada unidad en x que se avanza a través de esa línea tangente, se avanzan F(x0 , y0 ) unidades en y. Lo anterior quiere decir que por H unidades que se avanzan en x en la línea tangente se avanzarán (H)F(x0 , y0 ) unidades en y; y si la altura inicial es y0, la nueva altura sera: y1a = y0 + H F(x0 , y0 ) gráficamente:

Se usa y1a para denotar que esta altura es una aproximación a la altura real y1r de la solución, es decir solamente se esta calculando una aproximación a la curva solución. El error cometido en cada paso viene dado por e, en donde: e = y1r - y1a En la gráfica anterior, es evidente que entre más pequeño sea H más pequeño será e.

Si ahora trasladamos el problema de valor inicial al punto (x1 , y1a ) tenemos que se debe resolver el problema de valor inicial:

dy ---- = F(x, y) sujeto a y(x1 ) = y1a dx

para calcular la aproximación de la solución en el punto 2, es decir la coordenada (x2 , y2a ), gráficamente la situación es la siguiente:

O sea, la solución viene dada por: x2 = x1 + H y y2a = y1a + (H) F(x1 , y1a ) gráficamente (la interpretación es la misma):

Si continuamos haciendo el mismo tipo de razonamiento para i puntos, obtenemos la conocida fórmula de Euler: yia = y(i-1)a + H F(x(i-1) , y(i-1)a ) en donde: yia punto y(i-1)a H F(x(i-1) , y(i-1)a ) (x(i-1) , y(i-1)a )

=

Altura aproximada de la función solución desconocida en el siguiente

= = = =

Altura aproximada de la función solución desconocida en el punto actual Tamaño de paso Derivada de la función solución desconocida evaluada en el punto actual Coordenada del punto actual.

Inicialmente la coordenada del punto actual es (x0 , y0 ) —o sea la condición inicial. Si realizamos sistemáticamente la aplicación de la fórmula de Euler al problema de valor inicial planteado, lo que obtenemos es una gráfica aproximada —o bien las coordenadas aproximadas— de la función solución desconocida, es decir, resolvemos numéricamente la ecuación diferencial, gráficamente:

En la gráfica anterior, es evidente que si disminuimos el tamaño de paso (H) obtendremos una aproximación más precisa a la gráfica de la solución, pero a la vez requeriremos aplicar una mayor cantidad de veces la fórmula de Euler para obtener la solución en el mismo dominio o

intervalo ( x0 - xb ), veamos gráficamente el efecto de la reducción del tamaño de paso H en la aproximación de la solución:

Error computacional del método de Euler hacia adelante

Error de truncamiento Los errores de truncamiento son aquellos que resultan al usar una aproximación en lugar de un procedimiento matemático exacto (por ejemplo usar el método de Euler en vez de un método de integración analítica, para resolver una ecuación diferencial). El método de Euler hacia adelante debe utilizarse cuidadosamente para minimizar los errores de truncamiento. El error de truncamiento del método de Euler en un sólo intervalo [xi, xi+1 ] —llamado error local— es proporcional a H2, mientras que su error global —el correspondiente a considerar la totalidad de errores locales— al cubrir el dominio [x0 - xb ], es proporcional a H. Se deja como ejercicio al lector, usando el cálculo del error local absoluto presentado para el problema del Ejemplo 1, demostrar que estas aseveraciones son razonablemente correctas. El error de truncamiento es independiente del error computacional por redondeo. Consideremos la ecuación diferencial ordinaria de primer orden dy ---- = f(x,y), y(a) = c dx

El método de Euler para esta ecuación viene dado por: y(a + H ) = c + H f(a,c)

Sin embargo, suponiendo que f(x,c) sea diferenciable, el verdadero valor de y( a + H ) viene dado por la serie de Taylor: H2 y(a + H ) = y(c) + H f(a,c) + --- f '(a, c) + ... 2! comparando ambas ecuaciones es evidente que la diferencia entre ellas —el error cometido con el método de Euler— es H2 H3 e = --- f '(a,c) + ---- f '' (a,c) + ... 2! 3! Los valores f ' (a,c), f '' (a,c), ... se pueden encontrar diferenciando sucesivamente f(x,y) evaluada en el punto (a,c). Al usar el método de Euler, el error cometido, e, ocurre debido al hecho de que la serie se ha cortado —o truncado— en el segundo término —lo que explica su nombre de error de truncamiento—. Para valores de H suficientemente pequeños, se espera que el error fuera un número muy próximo al primer término —o sea a H2 / 2— ya que el resto de los términos aparecen divididos entre valores que crecen rápidamente (su tasa de crecimiento es factorial), siendo entonces cada vez más pequeños, por lo que se consideran despreciables. Pero he aquí una preguntita: — ¿Realmente lo son?—. Consideremos el uso del Teorema de Taylor con residuo. Este teorema establece que si f(x,y) es al menos doblemente diferenciable —o sea, si existe f '' (x,y)— entonces H2 y(a + H ) = y(c) + H f(a,c) + --- f '(r, c) 2! en donde r está entre a y (a + H)— o sea, a ≤ r ≤ H. El último término representa el residuo de la serie después de los primeros dos términos, superando el problema que tiene el Teorema de Taylor sin residuo de requerir todas las derivadas de f(x,y); por lo tanto el error pasa a ser:

H2 e = --- f '(r,c) 2! donde r esta entre los límites previamente especificados. Nótese que el error, e, es proporcional a H2, o en otras palabras, e es de orden H2 , simbólicamente e = O(H2 ), por lo que es evidente que reducir el tamaño de paso, H, reducirá considerablemente el error, e. Debido a que f '(r,c) puede ser positivo o negativo, si denotamos por M una constante positiva tal que |f ' (r,c)| < M, para valores de r entre a y (a + H ), entoces MH2 | E | < -----2 representando el término del lado derecho una cota superior para el error —un límite superior por encima del cual no estará nunca la magnitud del error absoluto (denotado por |E | ). Sin embargo este es el error cometido en un solo paso de "integración"o error local de truncamiento. Pero para obtener la solución completa requerimos n pasos de "integración" ¿Cuál sera el error acumulado en todos los pasos de "integración"? ¿De que orden sera la magnitud de éste error global de truncamiento? Para abordar estas dos preguntas consideremos lo siguiente: Formulemos la pregunta de esta manera: si tenemos y(a) y buscamos y(a + H) ¿Cuál sería el error acumulado si procedemos en n pasos al calculo de y(b), donde b = a + nH? Es posible profundizar el análisis del error para considerar el error acumulado, pudiendo demostrarse que el error máximo acumulado —o error global de truncamiento— es n veces el error absoluto dado antes, posiblemente con un valor diferente de M, digamos K; es decir, para n pasos de "integración" el error es KnH2 |En | < ------2 y dado que n = (b- a)/H K (b - a) H |En | < -------------2 lo cual permite ver que el error globlal de truncamiento del método de Euler es el del orden H, es decir en = O(H).

Error por inestabilidad

Un segundo tipo de error lo constituye la posible inestabilidad que aparece cuando la constante del tiempo, simbolizada por τ, es negativa. En este caso la solución numérica es inadecuada, a menos de que el intervalo de tiempo H sea suficientemente pequeño. Por ejemplo, una ecuación característica con solución decreciente —los valores de la variable independiente son cada vez menores— es dy ---- = -αy dx

s.a: y(0) = y0 > 0.

Su solución exacta es y = y0 e-αt El método de Euler para este problema es yn+1 = (1 - αH)yn si αH < 1, la solución numérica es decreciente y positiva; pero si αH > 2, la magnitud de la solución aumenta en cada paso y la solución oscila. Esto se conoce como inestabilidad [1 ] y plantea serios problemas de computación numérica al resolver cierto tipo de problemas; a este tipo de problemas se les conoce como ecuaciones rígidas[3 ] y se presentan frecuentemente al intentar analizar numéricamente problemas tales como la determinación de sistemas de flujo reactivo —por ejemplo al intentar evaluar la concentración del CO2 emitido por una chimenea, en diferentes puntos del tiempo y espacio. Generalmente αH < 1 para valores de H pequeños, por lo que es recomendable usar valores de H pequeños al utilizar el método de Euler hacia adelante. También se recomienda que H ≤ τ , en donde τ = α -1 es la constante de tiempo.

Error de redondeo

Los errores de redondeo se deben a que las computadoras solo guardan un número finito de cifras significativas durante un cálculo. Básicamente las computadoras redondean las cifras a fin de representarlas en memoria de alguna de estas dos formas: a) simplemente cortan los digitos sobrantes, o b) utilizan reglas de redondeo comúnes —lo cual agrega costo computacional adicional al problema—. Como quiera que lo hagan, este proceso se justifica bajo la suposición de que el número de cifras significativas en la mayor parte de las computadoras es mucho mayor que el error de redondeo dado por un corte —o redondeo verdadero— usualmente insignificante. Sin embargo, esta aseveración es discutible cuando se usan valores de H pequeños, que implican

un mayor número de cálculos que los requeridos para cubrir un dominio de solución dado, si se usaran valores de H mayores.

Error numérico total Si bien efectivamente, el reducir el tamaño de paso produce una solución más aproximada (ya que disminuye el error por truncamiento y por inestabilidad), también incrementa el número de operaciones que la computadora realiza para calcular la solución en un intervalo o dominio dado (incrementando consecuentemente el error de redondeo). El lector debe estar consciente de que existe un balance entre el tamaño mínimo de paso a utilizar (para minimizar el error de truncamiento) y el número máximo de operaciones a realizar (para minimizar el error de redondeo), el valor de H óptimo (aquel que minimiza el error numérico total) se encuentra entre ambos extremos. El usuario deberá ser cuidadoso al elegir el valor de H, o al juzgar la aplicabilidad del método de Euler a un cierto problema — existen métodos numéricos más avanzados basados en algoritmos de tamaño de paso variable y control de la cota máxima del error en cada paso, por ejemplo el método de Gear con control de paso y de error automático; sin embargo estos métodos avanzados se fundamentan primordialmente en el método de Euler. La siguiente gráfica ilustra la variación del error numérico total versus tamaño de paso (H); se ilustra tambien la variación de los errores de truncamiento y redondeo versus la misma variable independiente. Se muestra la localización del tamaño de paso óptimo, en donde los beneficios de su reducción en términos de error global de truncamiento son compensados por los perjuicios del excesivo número de cálculos computacionales en términos del error de redondeo.

Solución por computadora

El algoritmo del método de Euler ha sido repetidamente presentado en otras fuentes —ver por ejemplo[3 ]. Se presenta el código fuente de un programa en Pascal que implementa el método de Euler tal como ha sido explicado aquí: Program Euler; {Autor: Ing. Luis E. Castro-Solis. FIC UAdeC. Mexico 1996} Var x, y, x0, x1, y0, H, PI :real; NP, NC, i, j :longint;

NC = numero de pasos de calculo. }

Function F(x,y:real):real; Begin { programe aqui la EDO en la forma dy/dx = F(x,y) }

{ F:= 2*x*y

Diccionario de datos: End; x = valor actual de la var indep. y = valor actual de la var dep. x0 = valor inicial de la var indep. = inicio del dominio de soln. x1 = extremo final del dominio de soln. y0 = valor inicial de la var dep. H = tamaño de paso. PI = intervalo de impresion. NP = numero de pasos de impresion.

Begin {entrada de datos} Writeln; Writeln('Solucion de EDO''''s (Metodo de Euler).'); Writeln; Writeln('Datos...'); Writeln; Writeln('Dominio de la solución:'); Write('x0 = '); readln(x0); Write('Xf = '); readln(x1);

Writeln; Write('Condición inicial, y0 = '); readln(y0); Writeln; Write('Tamaño de paso, H = '); readln(H); Writeln; Write('Intervalo de salida tabulada = '); Readln(PI); {acondicionamiento de la salida impresa} NP := trunc( (x1 - x0) / PI); NC := trunc( PI / H ); {aplicacion del metodo de Euler} x

y := y0; Writeln; Writeln('Solución...'); Writeln; Writeln('-- x --':8,'- y --':8); Writeln(x:8:2, y:8:2); For i:=1 to NP do Begin For j:=1 to NC do begin y:= y + H*F(x,y); x:= x + H; end; writeln(x:8:2, y:8:2); end; Writeln; Writeln('ok.'); End.

:= x0;

El sofware anterior fué probado con éxito en la solución de los ejemplos resueltos presentados en la siguiente sección. El código anterior corré tal como esta bajo el compilador TURBOPascal v. 6.0 de la compañía Borland.

Ejemplos resueltos

Ejemplo 1

Resuelva y' = -20y + 7exp(-0.5t), y(0) = 5, por medio del método de Euler hacia adelante con h = 0.01, 0.001 y 0.0001 para 0 < t ≤ 0.09.

Solución Usando el código fuente presentado arriba en donde la función F se modificó de la siguiente manera para adaptarla a la ecuación diferencial del problema: Function F(x,y:real):real; Begin F:= -20*y + 7*exp(-0.5*x); End;

La salida impresa por el programa es la siguiente para h = 0.01 Solucion de EDO''s (Metodo de Euler). Datos... Dominio de la solucion: x0 = 0.00 Xf = 0.10 Condicion inicial, y0 = 5.00 Tamano de paso, H = 0.01 Intervalo de salida tabulada = 0.01

Solucion... -- x -- -- y -0.00 5.00 0.01 4.07 0.02 3.33 0.03 2.73 0.04 2.25 0.05 1.87 0.06 1.56 0.07 1.32 0.08 1.12 0.09 0.97 0.10 0.84 ok

para h = 0.001: Solucion de EDO''s (Metodo de Euler). Datos... Dominio de la solucion: x0 = 0.00 Xf = 0.10

Condicion inicial, y0 = 5.00 Tamano de paso, H = 0.00

Intervalo de salida tabulada = 0.01

0.03 0.04 0.04 0.05 0.06 0.07 0.08 0.09

Solucion... -- x -- -- y -0.00 5.00 0.01 4.23 0.02 3.58

3.04 2.60 2.22 1.91 1.65 1.43 1.25 1.10

ok

para h = 0.0001: Solucion de EDO''s (Metodo de Euler).

Solucion...

Datos...

-- x -- -- y -0.00 5.00 0.01 4.16 0.02 3.48 0.03 2.91 0.04 2.45 0.05 2.07 0.06 1.76 0.07 1.51 0.08 1.30 0.09 1.12 0.10 0.98

Dominio de la solucion: x0 = 0.00 Xf = 0.10 Condicion inicial, y0 = 5.00 Tamano de paso, H = 0.00 Intervalo de salida tabulada = 0.01

ok

Comentarios: Acerca del error numérico y la influencia del tamaño de paso. Consideremos nuevamente el problema de valor inicial: y' = -20y + 7exp(-0.5t)

s.a:

y(0) = 5

La expresión anterior es una EDO lineal de primer orden; la solución analítica de éste problema de valor inicial viene dada por: y = 5exp(-20t) + (7/19.5)[ exp(-0.5t) - exp(-20t) ] Se presenta abajo una gráfica de la solución analítica —"exacta"— y las soluciones numéricas tabuladas aproximadas con el programa Euler; Es posible apreciar la

diferencia entre la solución exacta y las soluciones analíticas computadas con el programa Euler a diferentes valores del tamaño de paso: 5,00 4,50 4,00 3,50 3,00 2,50 2,00 1,50 1,00 0,50 0,00 0,00

0,02

exacta

0,04

h = 0.01

0,06

h = 0.001

0,08

0,10

h = 0.0001

La exactitud del método de Euler hacia adelante aumenta al disminuir el tamaño de paso h. En efecto, las magnitudes de los errores son aproximadamente proporcionales a h. La siguiente tabla sumariza los resultados numéricos obtenidos contrastándolos contra la solución analítica exacta. Se incluye también el cálculo del error relativo absoluto asociado. Note como disminuye el error al disminuir el tamaño de paso.

Comparación de sol'ns calculadas con el programa Euler y la sol. analítica.

x

exact h = error a 0.01 0,00 5,00 5,00 0,000 0 0,01 4,16 4,07 0,020 9 0,02

3,47

3,33 0,039 3

0,03

2,90

2,73 0,058

Soluciones numéricas h= error h = error h = 0.001 0.0001 0.000001 5,00 0,000 5,00 0,0000 5,00 0 4,23 4,16 4,16 0,017 0,0007 6 3,58 3,48 3,47 0,032 0,0039 8 3,04 2,91 2,90

error 0,000 0 0,000 7 0,001 0 0,000

8 0,04

2,44

2,25 0,076 8

0,048 0 2,22 0,089 1

0,0032

0,05

2,06

1,87 0,091 1

1,91 0,071 7

2,07 0,0061

0,06

1,75

1,56 0,106 6

1,65 0,055 1

1,76 0,0079

0,07

1,49

0,08

1,28

0,09

1,11

1,43 0,041 0 1,25 0,024 9 1,10 0,009 3

0,1

0,97

1,32 0,114 7 1,12 0,126 3 0,97 0,126 4 0,84 0,133 6

1,51 0,0127 1,30 0,0141 1,12 0,0087 0,98 0,0108

2,45 0,0052

2 2,44 0,001 1 2,06 0,001 2 1,75 0,002 2 1,49 0,000 7 1,28 0,001 5 1,11 0,000 3 0,97 0,000 5

En la tabla también se presenta una solución numérica calculada con un tamaño de paso muy pequeño (H = 0.000001). Sin embargo, debe hacerse notar que una reducción mayor de H, sin el uso de números de computadora de doble precisión (tipo Double en TURBOPascal) tiene sus desventajas, puesto que aumenta el error numérico debido al redondeo.

Ejemplo 2

La tasa de enfriamiento de un cuerpo se puede expresar como: dT ---- = -k(T - Ta) dt en donde T es la temperatura del cuerpo (ºC), Ta es la temperatura del medio que rodea al cuerpo (ºC) y k es una constante de proporcionalidad (min-1 ). Por lo tanto, esta ecuación especifica que el enfriamiento de un cuerpo es proporcional a la diferencia de temperaturas entre el cuerpo y el medio que lo rodea. Si se calienta una bola de metal a 90 ºC y se sumerge en el agua que se mantiene a una

temperatura constante Ta = 20 ºC, empléese un método numérico para calcular el tiempo que le toma a la bola enfriarse a 30 ºC, sabiendo que k = 0.1 min-1 Solución La ecuación diferencial a resolver es:

dT ---- = 0.1(T - 20) s.a: T(0) = 90 dt adaptando el software al problema: Function F(x,y:real):real; Begin F:= 0.1*(y-20); End; la salida computada es: Solucion de EDO''s (Metodo de Euler). Datos... Dominio de la solucion: x0 = 0.00 Xf = 20.00 Condicion inicial, y0 = 90.00 Tamano de paso, H = 0.01 Intervalo de salida tabulada = 1.00

Solucion... -- x -- -- y -0.00 90.00 1.00 83.34 2.00 77.31 3.00 71.85 4.00 66.91 5.00 62.45 6.00 58.41 7.00 54.75 8.00 51.44 9.00 48.45 10.00 45.74 11.00 43.29 12.00 41.07 13.00 39.06 14.00 37.25 15.00 35.61 16.00 34.12 17.00 32.78 18.00 31.56 19.00 30.46

20.00

29.46

ok

Por lo que el tiempo solicitado es aproximadamente 19 minutos para que la bola se enfríe a 30 grados.

Ejemplo 3 Se disuelven inicialmente 50 libras (lb) de sal en un gran tanque que contiene 300 galones (gal) de agua. Se bombea salmuera al tanque a razón de 3 gpm; y luego la solución adecuadamente mezclada se bombea fuera del tanque también a razón de 3 gpm. Si la concentración de la solución que entra es de 2 lb/gal, determine la cantidad de sal que hay en el tanque en un instante cualquiera a) ¿Cuánta sal hay después de 50 min? b) ¿Cuanta después de un tiempo largo? Solución. Considerado al tanque completamente mezclado, homogéneo e isótropo, la cantidad de sal (lb) que hay en el tanque en un instante cualquiera viene dada por A(t ). La rapidez neta con que A(t ) cambia esta dada por — aplicando la Ley de Conservación de Masa— por el balance másico: dA ---- = (rapidez de entrada) - (rapidez de salida) = R1 - R2 dt La rapidez con que la sal entra al tanque es: R1 = ( 3 gal/min ) (2 lb/gal) = 6 lb/min La rapidez con la que la sal sale del tanque es: R2 = (3 gal/min) ( A(t )/300 lb/gal) = A(t )/100 lb/min En consecuencia la ecuación del balance de masa se transforma en: dA ---- = 6 - A(t )/100 dt la que se resolverá sujeta a la condición inicial A(0) = 50 lb Usando el software suministrado adaptado al problema:

Function F(x,y:real):real; Begin F:= 6 - y/100; End; la solución computada es la siguiente: Solucion de EDO''s (Metodo de Euler). Solucion... Datos... Dominio de la solucion: x0 = 0.00 Xf = 400.00 Condicion inicial, y0 = 50.00 Tamano de paso, H = 0.001 Intervalo de salida tabulada = 25.00

-- x -- -- y -0.00 50.00 25.00 171.66 50.00 266.41 75.00 340.20 100.00 397.67 125.00 442.42 150.00 477.28 175.00 504.43 200.00 525.57 225.00 542.03 250.00 554.85 275.00 564.84 300.00 572.62 325.00 578.67 350.00 583.39 375.00 587.07 400.00 589.93 ok

a) Para t = 50 se tiene que A = 266.41 lb. b) La solución analítica exacta del problema viene dada por: A(t ) = 600 - 550 exp(-t /100) de donde se puede ver que cuando t → ∞ , A(t ) → 600. Si se resuelve el problema con el programa presentado para un dominio de, digamos, 0 a 1000 se observará esta misma propiedad de la solución. Interpretando: el tanque alcanza un estado estable, es decir si no se cambia el régimen de alimentación la cantidad de sal que existirá en el tanque en forma estable será de 600 lbs, y a la larga el 2do término de la solución no tiene importancia.

Ejemplo 4

La ecuación diferencial dx --- = k x (n+1 - x ), k>0 dt

sujeta a la condición inicial x(0) = 1

proporciona un modelo razonable para describir la propagación de una epidema ocasionada al introducirse un individuo infectado en una población estática. La solución x( t ) representa el número de individuos infectados con la enfermedad en un instante cualquiera. Los sociólogos han utilizado este modelo para estudiar la difusión de información y el impacto de la publicidad en centros de población [4] Supóngase que un estudiante portador de un virus de gripe o influenza regresa a un campus universitario aislado que contiene 1000 estudiantes. Si se supone que la rapidez con que el virus se propaga es proporcional no solo al número x de estudiantes contagiados, sino también, al número de alumnos no contagiados, determinar el número de estudiantes contagiados después de 6 días, si además se observa que después de 4 días el número de estudiantes infectados es de x(4) = 50.

Solución: Suponiendo que el Campus se mantiene cerrado (nadie sale de él) durante el transcurso de la enfermedad, queremos resolver el problema de valor inicial: dx ----- = k x (1001 - x ) s.a: dt

x(0) = 1;

Usando el software suministrado adaptado al problema: Function F(x,y:real):real; Begin F:= k*y*(1001 - y) End; El valor de k (la tasa de propagación de la enfermedad) deberá encontrarse primeramente, calibrando el modelo para diferentes valor de k, hasta que la predicción para el cuarto día sea de 50 estudiantes infectados (valor observado).

Después de intentar varias veces para diferentes valores de k, se encontró que un valor de k = 0.0009906 permite predecir 50 estudiantes infectados al 4º día, como se muestra a continuación: Solucion de EDO''s (Metodo de Euler).

Intervalo de salida tabulada = 1.00000000

Datos... Solucion... Dominio de la solucion: x0 = 0.00 Xf = 4.00 Condicion inicial, y0 = 1.00 Tamano de paso, H = 0.00010000

-- x -- -- y -0.00 1.00 1.00 2.69 2.00 7.21 3.00 19.17 4.00 49.99 ok

Por lo tanto, se usará ese valor para predecir la propagación de la infección en el dominio solicitado:

Solucion de EDO''s (Metodo de Euler). Datos... Dominio de la solucion: x0 = 0.00 Xf = 11.00 Condicion inicial, y0 = 1.00 Tamano de paso, H = 0.00010000 Intervalo de salida tabulada = 1.00000000

Solucion... -- x -- -- y -0.00 1.00 1.00 2.69 2.00 7.21 3.00 19.17 4.00 49.99 5.00 124.12 6.00 276.19 7.00 506.79 8.00 734.54 9.00 881.68 10.00 952.53 11.00 981.83 ok

Se observa que el 6º día el número de estudiantes contagiados es de 276; eventualmente todos los estudiantes estarán contagiados, por lo que el valor de x =

1000 es asíntota de la solución x( t); esto se puede apreciar mejor graficando los resultados

x (No)

1000

500

0 0

5

10

15

t (dia)

Se observa la típica curva en "S" del crecimiento logístico con un periodo de aceleración exponencial y un periodo de desaceleración que tiende asintóticamente al valor 1000.

Método de Euler hacia adelante para sistemas de ecuaciones diferenciales

El método de Euler hacia adelante es aplicable a un conjunto de EDO's de primer orden. Consideremos un conjunto de EDO's de primer orden dado por: y ' = f( y, z, t) s.a: z' = g( y, z, t) s.a:

y(0) = y0 z(0) = z0

El método de Euler hacia adelante para el sistema anterior se escribe como: yn+1 = yn + Hy'n = yn + H f(yn, zn, tn) zn+1 = zn + Hz'n = zn + Hg(yn, zn, tn) En conclusión, una ecuación diferencial ordinaria de orden superior, digamos n, se puede descomponer en un sistema simultáneo de n ecuaciones diferenciales de primer orden.

Referencias

[1] pp.

Nakamura, S. Applied Numerical Methods with Software, Prentice Hall Inc. 1Ed. 570 Englewood Cliffs. 1991.

[2] Chapra, S.C. y Canale, R.P. Numerical Methods for Engineers with Personal Computer Applications, Mc. Graw-Hill Book Co., 1 Ed. 641 pp. New York. 1987. [3]

Burden, R. y Faires, D. Numerical Analysis, PWS, 3rd Ed. 732 pp. Boston. 1985.

[4] Zill, D.G. A First Course in Differential Equations with Applications, PWS, 3rd ed. 516 pp. Belmont, CA. 1986