1 Método de Diferencias Finitas para Ecuaciones en Derivadas Parciales Ecuaciones Elípticas y Parabólicas

2 A.uxx+B.uxy+C.uyy+D.ux+E.uy+F.u=GIntroducción Una ecuación diferencial en la que aparecen dos o más variables independientes se llama ecuación en derivadas parciales (EDP). EDP de orden 2, lineal y de coeficientes constantes: A.uxx+B.uxy+C.uyy+D.ux+E.uy+F.u=G Clasificación: Si B2 – 4.A.C < 0 , la ecuación se llama elíptica Si B2 – 4.A.C = 0 , la ecuación se llama parabólica Si B2 – 4.A.C > 0 , la ecuación se llama hiperbólica

3 EDP Elípticas

4 Solución Numérica de la Ecuación de LaplaceDominio de la solución

5 Diferencias Finitas i, j+1 i-1, j i, j i+1, j i, j-1

6 para una malla de cuadrados resultaAproximación por diferencias finitas para la solución de la ecuación de Laplace para una malla de cuadrados resulta Esta ecuación puede representarse en forma gráfica mediante el siguiente «Stencil»:

7 Ejemplo: cálculo de las líneas de corriente del escurrimiento fluido

8

9

10

11 SUCESIVAS POSICIONES DEL OPERADOR

12 Sistema de Ecuaciones ResultanteFormulación Matricial del Sistema denominando u=Y

13 Sistema de Ecuaciones Resultante

14 Resolución del Sistema en MATLAB >> b=[-50;-60; -10; 0; -30; -45] b = -50 -60 -10 -30 -45 >> psi=inv(A)*b psi =

15 Solución del Sistema de Ecuaciones

16 Valores de la función de corriente Y en la grilla discreta

17

18 EDP parabólicas Las ecuaciones parabólicas se emplean para caracterizar problemas dependientes del tiempo y el espacio Un ejemplo: ECUACIÓN DE CONDUCCIÓN DE CALOR

19 Ecuación de Conducción del CalorSe puede usar la conservación de calor para desarrollar un balance de energía en un elemento diferencial de una barra larga y delgada aislada, considerando la cantidad de calor que se almacena en un periodo de tiempo Δt Se llega a: en la que a es el coeficiente de difusividad térmica, que depende del material de la barra Caliente Frio L Para unas dadas condiciones iniciales y de contorno la solución de esta EDP parabólica permite conocer la temperatura en cualquier posición de la barra y para cualquier instante : u(x,t) con 0 < x < L ; t > 0

20 Condiciones de ContornoCondición Inicial Temperatura de la barra para t=0 u(x, 0) = f(x) (constante o función de la posición) Temperatura de la barra en los extremos x=0 ; x=L u(0, t) =T u(L, t) = TL Condiciones de Contorno

21 Solución Numérica de la Ecuación del CalorEncontraremos la solución aproximada de Para algunos puntos del dominio -algunas secciones de la barra Discretización del dominio -algunos instantes t x

22 Método de Diferencias finitasDiscretización del dominio (Grilla Discreta) Condiciones de Contorno e Iniciales en el Dominio Discretizado Reemplazo de las derivadas parciales de la EDP por sus aproximaciones numéricas. Se obtiene una Ecuación en Diferencias (ED) EDP ED Aplicación de la ED a los puntos de la Grilla Discreta

23 Convergencia Consistencia y EstabilidadEDP F(x,y,u)=0 Solución: ED Gi,j(h,k,u)=0, para cada (i,j) Convergencia Consistencia Estabilidad: la diferencia entre la solución numérica y la solución exacta tiende a cero a medida que avanza el cálculo con una cantidad de pasos tendiente a infinito

24 Discretización del DominioPara elegir los puntos en los cuales calcularemos la solución aproximada de la EDP -Particionamos el espacio secciones de la barra separadas una distancia h=Dx xi = i . h , i=0, 1, …m Elegimos una partición para el tiempo k=Dt tj = j . k , j=0, 1, …n h t j+1 k t j u i,j t j - 1 x x x i - 1 i i+1

25 Aproximaciones Numéricas de las derivadas parciales de la EDPDerivada parcial primera Hacia adelante Hacia atrás

26 Aproximaciones Numéricas de las derivadas parciales de la EDPDerivada parcial segunda Centrada

27 Ecuación del Calor. Método ExplícitoDerivada Primera : hacia delante Derivada Segunda: centrada Ecuación en diferencias Error O(k+h2) l parámetro de Courant

28 Ecuación del Calor. Método ExplícitoEcuación en diferencias Stencil

29 Estabilidad Numérica del Método ExplícitoEl método explícito se puede expresar en forma matricial: con w(j): conjunto de las aproximaciones u para el paso del tiempo j. w(j-1): conjunto de las aproximaciones u para el paso del tiempo j-1 A : matriz tridiagonal

30 Estabilidad Numérica del Método ExplícitoPara que el Método sea estable debe cumplirse que el radio espectral Esto es equivalente a que Como esto condiciona los valores de h y k que se adopten para la discretización (no puede utilizarse cualquier combinación de h y k).

31 Ecuación del Calor. Método explícito. EjemploHallar la temperatura para t = 0.3 s de una barra de 1m cuyos extremos se mantienen a 20ºC y a 40ºC. La temperatura inicial de la barra es de 100ºC y el coeficiente a = 0.1. Tomar Dx = 0.2m y Dt = 0.1 s. Justificar la aplicabili- dad del método explícito. Parámetro de Courant

32 Stencil

33 Grilla Discreta Cond. de Contorno Cond. Iniciales Valores calculados20º 40º 63.75º 89.06º 91.25º 72.81º 0.2s 20º 40º 70º 95º 96.25º 77.5º 0.1s 20º 40º 80º 100º 100º 85º 0s 20º 40º 100º 100º 100º 100º 0m 0.2m 0.4m 0.6m 0.8m 1.0m Cálculo de Temperaturas con el Stencil: Valores calculados T0.2m,0.1s=0.25* * *100=80º T0.6m,0.1s=0.25* * *100=100º T0.4m,0.1s=0.25* * *100=100º T0.8m,0.1s=0.25* * *40=85º

34 Método Explícito. Cálculo con Excel

35 Ecuación del Calor. Método ImplícitoDerivada Primera : hacia atrás Derivada Segunda: centrada Ecuación en diferencias Error O(k+h2) l parámetro de Courant

36 Ecuación del Calor. Método ImplícitoEcuación en diferencias Stencil Este esquema permite plantear la solución en el instante tj a partir de la solución en el instante tj-1

37 Estabilidad Numérica del Método ImplícitoLa aplicación sucesiva del stencil da como resultado un sistema de ecuaciones, que resulta ser tridiagonal. El método implícito se puede expresar en forma matricial: con w(j): conjunto de las aproximaciones u para el paso del tiempo j. w(j-1): conjunto de las aproximaciones u para el paso del tiempo j-1 A : matriz tridiagonal

38 Estabilidad Numérica del Método ImplícitoAnalizando la Matriz A correspondiente al Método Implícito, se observa que el radio espectral Para cualquier valor del parámetro de Courant, por lo que este método es incondicionalmente estable.

39 Ecuación del Calor. Método implícito. EjemploHallar la temperatura para t = 0.3 s de una barra de 1m cuyos extremos se mantienen a 20ºC y a 40ºC. La temperatura inicial de la barra es de 100ºC y el coeficiente a = 0.1. Tomar Dx = 0.2m y Dt = 0.1 s. Parámetro de Courant

40 Stencil

41 Grilla Discreta Cond. de Contorno Cond. Iniciales 0.3s 20º 69.05º88.85º 90.55º 76.54º 40º T0.2,0.3 T0.4,0.3 T0.6,0.3 T0.8,0.3 0.2s 20º 40º 76.37º 93.37º 94.48º 82.17º T0.2,0.2 T0.4,0.2 T0.6,0.2 T0.8,0.2 0.1s 20º 86,22º 97,34º 97,83º 89,64º 40º T0.2,0.1 T0.4,0.1 T0.6,0.1 T0.8,0.1 0s 20º 40º 100º 100º 100º 100º 0m 0.2m 0.4m 0.6m 0.8m 1.0m -0.25* *T0.2, *T0.4, *86.22=0 -0.25* *T0.2, *T0.4, *100=0 -0.25* *T0.2, *T0.4, *76.37=0 -0.25*T0.2, *T0.4, *T0.6,0.2-1*97.34=0 -0.25*T0.2, *T0.4, *T0.6,0.3-1*93.37=0 -0.25*T0.2, *T0.4, *T0.6,0.1-1*100=0 -0.25*T0.4, *T0.6, *T0.8,0.2-1*97.83=0 -0.25*T0.4, *T0.6, *T0.8,0.3-1*94.48=0 -0.25*T0.4, *T0.6, *T0.8,0.1-1*100=0 -0.25*T0.6, *T0.8, *40 - 1*89.64=0 -0.25*T0.6, *T0.8, *40 - 1*82.17=0 -0.25*T0.6, *T0.8, *40 - 1*100=0

42 Método de Crank-NicolsonIdea: obtener un método con error O(k2+h2) ¿Cómo? Se promedian las diferencias hacia delante en el j-ésimo paso en t y las diferencias hacia atrás en el (j+1)ésimo paso en t El método es incondicionalmente estable Stencil