• Categories
• Ecuaciones
• Derivado
• Ecuaciones diferenciales
• Análisis numérico
• Linealidad

Citation preview

Universidad Nacional Agraria La Molina Facultad de Ingeniería Agrícola Departamento de Recursos Hídricos

Análisis Numérico para Ingeniería (Aplicaciones con MATLAB)

Ecuaciones Diferenciales Ordinarias Problemas de Valor Inicial Jesús Abel Mejía Marcacuzco, Ph.D. Lima - Perú

ECUACIONES DIFERENCIALES ORDINARIAS

El comportamiento dinámico de los sistemas implica desplazamientos, velocidades y aceleraciones y sus derivadas respecto al tiempo de estas cantidades. Una ecuación en la que intervienen una o más derivadas ordinarias de la función incógnita se denomina ecuación diferencial ordinaria (EDO). El orden de la ecuación está determinado por el orden de la derivada más alta. Los problemas de EDO se clasifican en problemas de valor inicial y problemas de valor en la frontera, dependiendo de cómo se especifican las condiciones en los extremos del dominio. Todas las condiciones de un problema de valor inicial se especifican en el punto inicial. El problema se convierte en uno de valor en la frontera si las condiciones se extienden entre los puntos inicial y final. Las EDO en el dominio del tiempo son problemas de valor inicial, así que todas las condiciones se especifican en el momento inicial, como t = 0.

PROBLEMAS DE VALOR INICIAL EDO de Primer Orden El problema de valor inicial de una EDO de primer orden puede escribirse en la forma:

y' (t )  f  y, t ,

y0  y0

donde f(y,t) es función de y y t; con y(0) = y0 como condición inicial indispensable para la solución. La primera derivada de y se da en función de y y t, y lo que se quiere calcular es la función desconocida y integrando numéricamente f(y,t). Si f fuera independiente de y, el cálculo sería un problema de integración directa. Si f es una función lineal de y, por ejemplo: f = ay + b, con a y b son constantes o funciones de t, la ecuación es una EDO lineal. Si f es una función no lineal de y, la ecuación es una EDO no lineal. Es posible encontrar la solución analítica de algunas EDO, pero la mayor parte de ellos son no lineales y carecen de solución analítica y debemos usar métodos numéricos.

Ejemplo: Un cuerpo de masa M se deja caer desde un avión en t = 0. Si la velocidad inicial es cero y la caída vertical; deduzca una EDO para determinar la velocidad vertical del cuerpo, considerando la fuerza de arrastre aerodinámico igual a Faire = CV2. Solución: Aplicando la primera ley de Newton, con V como la velocidad vertical (positiva hacia abajo) y con g como la aceleración de la gravedad; se obtiene: dV M  gM  Faire dt

La EDO obtenida es:

La ecuación puede escribirse en términos de y que es la distancia recorrida por el cuerpo en su caída:

dV C 2 g V dt M

d2y C  2 dt M

2

 dy    g  dt 

Métodos de Euler hacia adelante

Los métodos de Euler son procedimientos sencillos para resolver EDO de primer orden y se pueden programar con gran facilidad. Los métodos de Euler incluyen tres versiones: hacia delante, modificado y hacia atrás y su exactitud no es muy alta. El método de Euler hacia adelante para y’ = f(y,t) se deduce rescribiendo la aproximación de diferencia hacia delante: y n 1  y n  y n' h

y n1  y n  hf  y n , tn   hyn'

Recursivamente la ecuación puede representarse como:

y1  y0  hy0'  y0  hf  y0 ,t0 

y 2  y1  hf  y1 ,t1 

y3  y 2  hf  y2 ,t 2  . y n  yn1  hf  y n1 , t n1 

Ejemplo:

Usar el método de Euler para integrar numéricamente la siguiente ecuación.

dy  2 x 3  12 x 2  20 x  8.5 dx

Comparar los resultados con la solución analítica y calcular el error relativo porcentual

y  0.5x 4  4x3 10x2  8.5x  1

x 0.000 0.500 1.000 1.500 2.000 2.500 3.000 3.500 4.000

y y y (h=0.5) (h=0.25) (h=0.125) Exacto 1.000 1.000 1.000 1.000 5.250 4.180 3.686 3.219 5.875 4.344 3.648 3.000 5.125 3.555 2.857 2.219 4.500 3.125 2.531 2.000 4.750 3.617 3.139 2.719 5.875 4.844 4.398 4.000 7.125 5.867 5.279 4.719 7.000 5.000 4.000 3.000

Error (h=0.5) 0.000 -63.107 -95.833 -130.986 -125.000 -74.713 -46.875 -50.993 -133.333

Error (h=0.25) 0.000 -29.854 -44.792 -60.211 -56.250 -33.046 -21.094 -24.338 -66.667

Error (h=0.125) 0.000 -14.502 -21.615 -28.785 -26.563 -15.445 -9.961 -11.879 -33.333

Ejemplo:

Usar el método de Euler para integrar numéricamente la siguiente ecuación.

dy 2 t  y dt

Comparar los resultados con la solución analítica y calcular el error relativo porcentual

y  et  t 2  2t  2

t 0.000 0.400 0.800 1.200 1.600 2.000 2.400 2.800 3.200

y y y Error Error Error (h=0.4) (h=0.2) (h=0.1) Exacto (h=0.4) (h=0.2) (h=0.1) 1.000 1.000 1.000 1.000 0.000 0.000 0.000 0.700 0.648 0.670 0.690 -1.496 6.043 2.925 0.538 0.512 0.553 0.591 8.917 13.265 6.449 0.569 0.630 0.686 0.739 23.038 14.689 7.173 0.830 1.026 1.093 1.158 28.329 11.426 5.602 1.349 1.714 1.791 1.865 27.654 8.075 3.973 2.144 2.705 2.788 2.869 25.267 5.725 2.826 3.229 4.005 4.093 4.179 22.736 4.172 2.065 4.612 5.617 5.709 5.799 20.467 3.134 1.554

Ejemplo

a) Resolver la siguiente EDO utilizando el método de Euler hacia adelante con h = 0.01 para 0 < t  0.02. Haga los cálculos de esta parte a mano. y '  f  y, t , y 0  5 donde f  y, t   20 y  7 exp  0.5t  b) Repita para 0 < t  0.1 con h = 0.01, 0.001 y 0.0001. Evalúe los errores de los tres cálculos comparándolos con la solución analítica. Solución

a) Los primeros dos incrementos de tiempo de los cálculos con h = 0.01 son:

t0  0, t1  0.01, t2  0.02,

y0  y0  5

y1  y0  hy0'  5  0.01 205  7 exp 0  4.07 y2  y1  hy1'  4.07  0.01 204.07  exp 0.01  3.326

b) La solución analítica se obtiene de resolver la ecuación: 20t dy 20t  0.5t 20t dy  0.5 t Multiplicando por e 20t e  20 ye  7 e e  20 y  7e dt dt  7 0.5t 90.5  20t  7 19.5t 20t d ye 20t 19.5t y e  e  ye  e C  7e 19.5  19.5  19 .5 dt





C = 90.5/19.5 calculado para t = 0; y =5 t 0.00 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.10

h=0.01 5.00 4.07000 3.32565 2.72982 2.25282 1.87087 1.56497 1.31990 1.12352 0.96607 0.83977

h=0.001 5.00 4.14924 3.45379 2.88524 2.42037 2.04023 1.72932 1.47496 1.26683 1.09646 0.95696

h=0.0001 5.00 4.15617 3.46513 2.89915 2.43554 2.05574 1.74454 1.48949 1.28041 1.10895 0.96830

Exacto 5.00000 4.15693 3.46638 2.90068 2.43721 2.05745 1.74622 1.49109 1.28191 1.11033 0.96956

% Error % Error h=0.01 h=0.001 0.000 0.000 2.091 0.185 4.060 0.363 5.890 0.532 7.566 0.691 9.069 0.837 10.380 0.968 11.481 1.082 12.356 1.176 12.993 1.249 13.386 1.300

% Error h=0.0001 0.000 0.018 0.036 0.053 0.069 0.083 0.096 0.108 0.117 0.124 0.130

La exactitud del método de Euler aumenta al reducirse el tamaño del incremento h.

Métodos de Euler modificado El método de Euler modificado es más exacto y más estable que el método hacia delante. El método de Euler modificado se deduce aplicando la regla trapezoidal a la solución de y’ = f(y,x): h yn 1  yn   f  yn 1 , tn 1   f  yn , tn  2 Si f es una función lineal de y, es fácil resolver la ecuación para f  y, t   ay  cost  yn+1 en forma cerrada. Por ejemplo si h Se obtiene para yn+1: yn 1  yn  ayn 1  cos tn 1   ayn  cos tn  2 Si f es una función no lineal de y, la ecuación se convierte en una función no lineal de yn+1, por lo que debemos aplicar un algoritmo para resolver la ecuación no lineal. Un método muy utilizado es el de sustituciones sucesivas, que se escribe así: y n 1 

1  ah / 2 h cos tn 1  cos tn  yn  1  ah / 2 1  ah / 2 2

Si f es una función no lineal de y, la ecuación se convierte en una función no lineal de yn+1, por lo que debemos aplicar un algoritmo para resolver la ecuación no lineal. Un método muy utilizado es el de sustituciones sucesivas, que se escribe así: k 

y n 1





h  y n  f  y nk11 , t n 1   f  y n , t n  2

k  Donde yn 1 es la k-ésima aproximación iterativa para yn 1 y

yn01

k  es una estimación inicial de yn 1 . La iteración termina cuando yn 1

es menor que

k 1 para para una tolerancia especificada. y n1

Si tomamos como estimación inicial yn, el primer paso de la iteración será idéntico al método de Euler hacia delante. Si sólo se emplean dos pasos de iteración, el método se convierte en el método de Runge-Kutta de segundo orden.

Ejemplo a) Para el siguiente problema de valor inicial determine y(0.1) por el método de Euler modificado con h = 0.1. Establezca como tolerancia de la convergencia 0.00001. y'   y1.5  1, y 0  10 b) Continúe el cálculo avanzando incrementos de tiempo hasta llegar a t = 0.5 c) Elabore un algoritmo en MATLAB que calcule la solución para 0 < t  1 por los métodos de Euler hacia delante y modificado y grafique los resultados.



Solución a) El esquema de Euler modificado se escribe:



yn 1  yn  h / 2   yn 1    yn   2



1.5

1.5





y su solución iterativa basada en sustituciones sucesivas es:

yn 1  yn  h / 2  yn 1 k 

k 1

  y  1.5

Donde k es el número de iteración.

1.5

n

2

La iteración para y1 comienza con la estimación inicial y continua así:

y10  y0  10,





y1  10  0.1 / 2  10  10  2  6.93772 1

1.5

1.5

 y    10  0.1 / 2 7.60517 y    10  0.1 / 2 7.47020 y    10  0.1 / 2 7.49799



y12   10  0.1 / 2  6.93772  10  2  7.60517 1.5

1.5

  2  7.49799  2  7.49229

3

1.5

 10  2  7.47020

4

1.5

 10

1

1

1.5

5

1

...



1.5

1.5

 10

1.5



y1  10  0.1 / 2  7.49326  10  2  7.49326 1.5

1.5

b) El cálculo hasta t = 1.0, se muestra a continuación: Y(t=0) 0.0 10 6.93772 7.60518 7.47020 7.49800 7.49229 7.49346 7.49322 7.49327 7.49326 7.49326

Y(t=0.1) 0.1 7.49326 5.54207 5.91532 5.84832 5.86051 5.85830 5.85870 5.85863 5.85864 5.85864 5.85864

t 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0

Y(t=0.2) 0.2 5.85864 4.54058 4.76584 4.72940 4.73535 4.73438 4.73454 4.73451 4.73452 4.73452 4.73452

y 10.0 7.49326 5.85864 4.73452 3.92990 3.33574 2.88599 2.53861 2.26584 2.04868 1.87382

Y(t=0.3) 0.3 4.73452 3.80434 3.94841 3.92714 3.93031 3.92983 3.92990 3.92989 3.92990 3.92990 3.92990

Y(t=0.4) 0.4 3.92990 3.25083 3.34730 3.33416 3.33596 3.33571 3.33575 3.33574 3.33574 3.33574 3.33574

Y(t=0.5) 0.5 3.33574 2.82650 2.89353 2.88502 2.88611 2.88597 2.88599 2.88599 2.88599 2.88599 2.88599

Y(t=0.6) 0.6 2.88599 2.49571 2.54371 2.53800 2.53868 2.53860 2.53861 2.53861 2.53861 2.53861 2.53861

Y(t=0.7) 0.7 2.53861 2.23413 2.26940 2.26543 2.26588 2.26583 2.26584 2.26584 2.26584 2.26584 2.26584

Y(t=0.8) 0.8 2.26584 2.02477 2.05124 2.04841 2.04871 2.04868 2.04869 2.04868 2.04868 2.04868 2.04868

Y(t=0.9) 0.9 2.04868 1.85545 1.87570 1.87362 1.87384 1.87382 1.87382 1.87382 1.87382 1.87382 1.87382

Y(t=1.0) 1.0 1.87382 1.71732 1.73304 1.73149 1.73165 1.73163 1.73163 1.73163 1.73163 1.73163 1.73163

c) El Archivo-M para la solución del problema en MATLAB y los resultados, se muestran a continuación. En el guión, la solución del método de Euler hacia adelante se denota con yf, mientras que la del método modificado es ym. La ecuación de Euler modificada se resuelve iterativamente por el método de sustitución sucesiva con el número de iteraciones limitado a 10 como máximo. Si el número de iteración rebasa 9, se imprime un mensaje. clear, clf, hold off %===== Euler hacia adelante yf(1)=10; t(1)=0; h=0.1; n=1; while t(n)