Document 7324901

Download Report

Transcript Document 7324901

Métodos de Diferencias Finitas
para
Ecuaciones en Derivadas Parciales
Ecuaciones en Derivadas Parciales

Introducción

Diferencias finitas

Convergencia y estabilidad

Ecuaciones hiperbólicas: ecuación de Ondas

Ecuaciones parabólicas: ecuación del Calor

Ecuaciones elípticas: ecuación de Laplace
Introducción


EDP de orden 2, lineales de coeficientes
constantes.
Auxx+Buxy+Cuyy+Dux+Euy+Fu=G
 Ecuación
de Ondas
utt - c2uxx = 0
 Ecuación
del Calor
ut - cuxx = 0,
 Ecuación
de Laplace
uxx + uyy = 0
Condiciones iniciales y de contorno
c>0
Diferencias finitas

Discretización:

Métodos explícitos
 Sencillos
 Inestables

Métodos implícitos
 Más
complejos
 Estables
EDP
EDF
h
yj+1
k
yj
u i,j
yj - 1
x i-1
x i x i+1
Diferencias primeras

Hacia adelante
u x ( xi , y j ) 

h

u i +1, j - u i , j
h
Error
u x ( xi , y j ) 

u( x i + h, y j ) - u( x i , y j )
u i +1, j - u i , j
h
h
- u xx ( x i + h, y j ), 0    1
2
Hacia atrás
u x ( xi , y j ) 
u( x i , y j ) - u( x i - h, y j )
h

u i , j - u i -1, j
h
Diferencias primeras (cont.)

Diferencias simétricas
u x ( xi , y j ) 

u( x i + h, y j ) - u( x i - h, y j )
2h

u i +1, j - u i -1, j
2h
Error
u( x + h, y) - u( x - h, y) h 2
u x ( x, y) 
u xxx (, y),
2h
6
 ]x - h, x + h[
Diferencias segundas

Diferencias simétricas
u xx ( x i , y j ) 

u x ( x i + h, y j ) - u x ( x i , y j )
h

u i +1, j - 2 u i , j + u i -1, j
h2
Error
u( x - h, y) - 2 u( x, y) + u( x + h, y) h
u xx ( x, y) 
- u xxxx (, y)
2
k
12
 ]x - h, x + h[
2
Convergencia y estabilidad
Solución:
~
u ( x, y)

EDP
F(x,y,u)=0

EDF
Gi,j(h,k,u)=0, para cada (i,j) u h ,k ( x i , y j )

Convergencia
u h ,k ( xi , y j )
h
0 u( xi , y j )
,k
G i , j ( h, k , ~
u)

Consistencia

Estabilidad:

Consistencia + Estabilidad
h
0 0
,k 
Control del error de redondeo
Convergencia
Ecuaciones hiperbólicas

Ecuación de
Ondas
utt = c²uxx , 0 < x < L, t > 0

Condiciones
iniciales
u(x, 0) = f(x)
ut(x, 0) = g(x)

Condiciones
de contorno
u(0,t) = l(t)
u(L,t) = r(t)

Ecuación en diferencias finitas
ui , j +1 - 2 ui , j + ui , j -1
u
- 2ui , j + ui -1, j
2 i +1, j
c
2
k
h2
Ec. de Ondas: Método explícito

Condiciones iniciales
ui,0 = fi

y
ui,1 - ui,-1 = 2kgi
Paso 1º
ui,1 = a2 (fi-1+fi+1)/2 + (1-a2)fi + kgi

Pasos siguientes
ui,j+1 = a2(ui+1,j + ui-1,j) +2(1 - a2)ui,j - ui,j-1

Convergencia
a1
Ecuación de ondas. Método explícito.
utt = c²uxx , 0 < x < L, t > 0
Ejemplo
c = 1, L=T=4,
u(x, 0) = 2-|x-2|
nx=4, nt=8,
ut(x, 0) = 0
u(0,t) = 0
u(L,t) = 0
Condición de convergencia : a  ck
h
 Instante t = 0:


1 0.5 1
 1
1
2
u0,0 = f(x0) = 2 - |x0 - 2| = 2 - |0 - 2| = 0 = f(x4)
u1,0 = f(x1) = 2 - |x1 - 2| = 2 - |1 - 2| = 1 = f(x3)
u2,0 = f(x2) = 2 - |x2 - 2| = 2 - |2 - 2| = 2

Instante t=1:
ui,1 = a2·(ui-1,0+ui+1,0)/2 + (1 - a2)·ui,0 + k·g(xi)
donde a2 = 1/4, 1 - a2 = 3/4:
u1,1 = (1/4)(u0,0 + u2,0)/2 + (3/4)u1,0 =
(1/4)(0 + 2)/2 + (3/4)1 = 1 = u3,1
u2,1 = (1/4)(u1,0 + u3,0)/2 + (3/4)u2,0 =
(1/4)(1 + 1)/2 + (3/4)2 = 7/4
Aplicando la fórmula genérica
ui,j+1 = a2·(ui-1,j + ui+1,j) + 2·(1 - a2)·ui,j - ui,j-1
con lo que, para t = 1 obtenemos:
u1,2 = (1/4)(u0,1 + u2,1) + (3/2)u1,1 - u1,0
= (1/4)(0 + 7/4) + (3/2)1 - 1 = 15/16
= u3,2
u2,2 = (1/4)(u1,1 + u3,1) + (3/2)u2,1 - u2,0
= (1/4)(1 + 1) + (3/2)(7/4) - 2 = 9/8
Procediendo análogamente
t=0
t = 0.5
t=1
t = 1.5
t=2
t = 2.5
t=3
t = 3.5
t=4
x=0
0
0
0
0
0
0
0
0
0
x=1
x=2
1.0000
2.0000
1.0000
1.7500
0.9375
1.1250
0.6875
0.4063
0.1953 -0.1719
-0.4375 -0.5664
-0.9932 -0.8965
-1.2764 -1.2749
-1.2401 -1.6541
x=3
1.0000
1.0000
0.9375
0.6875
0.1953
-0.4375
-0.9932
-1.2764
-1.2401
x=4
0
0
0
0
0
0
0
0
0
Ec. de Ondas: Método implícito
Idea
ui,j+1 - 2ui,j + ui,j-1 = a2 [(ui+1,j+1 - 2ui,j+1 + ui-1,j+1)

+ (ui+1,j-1 - 2ui,j-1 + ui-1,j-1)]/2

Pasos
(1+a2)ui,j+1 - a2(ui+1,j+1 + ui-1,j+1)/2 =
2ui,j + a2(ui+1,j-1 + ui-1,j-1)/2 - (1+a2)ui,j-1

Convergencia
para todo a
Algoritmo del método implícito

Truco ecuación implícita
- a2( ui-1,j-1 + ui-1,j+1)/4
+ (1 + a2)(ui,j-1 + ui,j+1)/2
- a2(ui+1,j-1 + ui+1,j+1)/4 = ui,j

Sistema
tridiagonal

Factorización LU
.
Aw = v,
v = (u1,j,u2,j,...,unx-1,j)'
ui,j+1 = wi - ui,j-1
Lz = v
Uw = z
Método implícito.

Resolución del sistema

Sustitución
 58

 -116
 0

-1
16
5
8
-1
16
0  x1   1 
   
7
-1
16  x1    4 
1
5  x 
8  1 
 
 u1, 2   x1   u1, 0   1.9184   1   0.9184 

   
 
   

 u2, 2    x1  -  u2, 0    3.1837  -  2    1.1837 
 u   x   u   1.9184   1   0.9184 
   

 3 , 2   1   3, 0  

Factorización LU
0
0 
 0.625


L   - 0.0625 0.61875
0 


0
0
.
0625
0
.
61
86


0 
 1 - 0 .1


U  0
1
- 0.10 
0

0
1


t=0
t = 0.5
t=1
t = 1.5
t=2
t = 2.5
t=3
t = 3.5
t=4
x=0
0
0
0
0
0
0
0
0
0
x=1
1.0000
1.0000
0.9184
0.6926
0.2912
-0.2449
-0.7996
-1.2231
-1.3966
x=2
2.0000
1.7500
1.1837
0.4824
-0.1699
-0.6647
-0.9953
-1.2214
-1.3981
x=3
1.0000
1.0000
0.9184
0.6926
0.2912
-0.2449
-0.7996
-1.2231
-1.3966
x=4
0
0
0
0
0
0
0
0
0
Ecuaciones parabólicas
Ecuación
ut = cuxx, 0 < x < L, t > 0
del Calor
 Condición
u(x, 0) = f(x)
inicial
 Condiciones
u(0, t) =T0 u(L, t) = TL
de contorno
 Ecuación en diferencias

ui , j + 1 - ui , j
ui +1, j - 2ui , j + ui -1, j
c
k
h2
Ec. del Calor: Método explícito

Condición inicial
ui,0 = f(xi)

Condiciones de contorno
u0,j = T0

unx,t = TL
para j>0
Pasos siguientes
ui,j+1 = a(ui+1,j+ui-1,j) +(1-2a)ui,j

Convergencia
a  1/2
Óptimo a = 1/6
Ecuación del Calor. Método explícito.
Ejemplo

Hallar la temperatura para t = 0.3 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 c = 0.1. Tomar Dx = 0.2 y
Dt = 0.1. Justificar la aplicabilidad del método
explícito.
ck 0.1  0.1 1 1
 a 2 
 
h
0.04
4 2
Ajuste de las condiciones iniciales y de
contorno:
u0,0 = 60, u1,0 = u2,0 = u3,0 = u4,0 = 100,
u5,0 = 70
 Instante t = 0.1
u1,1 = (u0,0 + u2,0)/4 + u1,0/2
= (60+100)/4 + 100/2 = 90
u2,1 = u3,1 = 100
u4,1 = (u3,0 + u5,0)/4 + u4,0/2
= (100+70)/4 + 100/2 = 92.5



Instante t = 0.2 :
u1,2 = 75
u3,2 = 98.125
u2,2 = 97.5
u4,2 = 81.25
Instante t = 0.3:
u1,3 = 66.875
u3,3 = 93.75
u2,3 = 92.0313
u4,3 = 75.1563
Ec. del Calor: Método implícito

Idea: Diferencias hacia atrás
ui , j - ui , j -1
ui +1, j - 2ui , j + ui -1, j
c
2
k
h

Pasos
(1+2a)ui,j - a(ui-1,j + ui+1,j) = ui,j-1

Convergencia
para todo a
Ecuación del Calor. Método implícito
Se verifica la condición de convergencia :
a = 1/4 < 1/2
 Diagonal principal: 1 + 2a = 3/2,
Diagonales contiguas -a = -1/4.
 Para t = 0.1: 3
 2 - 14
 u1,1   u1, 0 + a u0, 0 



  
u2, 0

 - 14 3 2 - 14
 u2,1  






3
1
1
u
u
- 4
3, 0
2
4


 3,1  

 u +a u 
3  u
1
4
,
1
5, 0 
4
2

   4, 0

Valores obtenidos por este método:
t = 0.1
t = 0.2
t = 0.3
x = 0.2
86.2237
76.3776
69.0598
x = 0. 4 x = 0.6
97.3423 97.8301
93.3707 94.4771
88.8487 90.5494
x = 0.8
89.6384
82.1718
76.5394
Método de Crank-Nicholson

Idea: media de diferencias centrales
ui,j+1 - ui,j = a [(ui+1,j+1 - 2ui,j+1 + ui-1,j+1) +
(ui+1,j - 2ui,j + ui-1,j)] /2

Pasos
2(1+a)ui,j+1 - a(ui+1,j+1 + ui-1,j+1) =
2(1-a)ui,j + a(ui+1,j + ui-1,j)

Convergencia
para todo a
Ecuación del Calor. Método de
Crank-Nicholson
 5 2 - 14



5
1
1
- 4
- 4

2

- 14 5 2 - 14 



5 
1
4
2 


Matriz del sistema:

Término independiente del primer paso:
100 
 60 + 100 




3 100  1 100 + 100 
+ 


2 100 4 100 + 100 




100 
 100 + 70 





Valores obtenidos por Crank-Nicholson:
t = 0.1
t = 0.2
t = 0.3
x = 0.2
x = 0.4
87.8683 98.6826
76.0999 95.1069
68.2003 90.2963
x = 0.6
x = 0.8
98.9578 90.8958
96.0470 82.0380
91.9748 76.0250
Ecuaciones elípticas

Ecuación de Laplace
uxx + uyy = 0,

0 < y <b
Condiciones de contorno
u(x,0),

0 < x < a,
u(x,b),
u(0,y),
u(a,y)
Discretización
u i+1, j - 2 u i , j + u i-1, j
h
2
+
u i , j+1 - 2u i , j + u i , j-1
k
2
0
Ecuación de Laplace

Ecuación en diferencias: a=k/h
a2(ui-1,j + ui+1,j) + ui,j-1 + ui,j+1 - 2(a2+1)ui,j = 0

Matriz del
sistema:
grande ,
dispersa

Caso h = k :
ui-1,j + ui+1,j + ui,j-1 + ui,j+1 = 4ui,j
Ec. de Laplace: Métodos iterativos

Método de Jacobi
 

u (i ,kj)  a 2 u (i -k1-,1j) + u (i +k1-,1j) + u (i ,kj--11) + u (i ,kj+-11)


2 + 2a2 
Método de Gauss-Seidel
 

u(i ,kj)  a 2 u(i -k1) , j + u(i +k1-,1j) + u(i ,kj)-1 + u(i ,kj-+11)


2 + 2a 2 
Criterio de parada
max u (i ,kj) - u (i ,kj-1)  
i, j
Método de Sobrerrelajación
Idea:
ponderar el desplazamiento de Gauss-Seidel
 Pasos

 

u (i ,kj)  a 2 u(i -k1) , j + u(i +k1-,1j) + u(i ,kj)-1 + u(i ,kj-+11)
u

(k)
i, j
 (1 - w ) u
( k - 1)
i, j

(k)

+ wui , j
Si w = 1 coincide con Gauss-Seidel
2 + 2a 2 
Ecuación de Laplace.
Ejemplo
uxx+ uyy=0, 0 < x < 1
0 < y < 1,
n=4
m=4,
u(x, 0) = 0
u (x, 1) = 100x
u(0, y) = 0
u(1, y) = 100y

 = 0.01

Ajuste de las condiciones de contorno:
ui(,0j )  25
i  1,...,m-1
j 1,...,n-1
Método de Jacobi.
Iteraciones: 8
 Operaciones en coma flotante: 1142

x = 0.25
x = 0.5
x = 0.75
y = 0.25 y = 0. 5 y = 0.75
6.2561 12.5031 18.7500
12.5031 25.0000 37.4969
18.7500 37.4969 56.2439
Método de Gauss-Seidel.
Iteraciones: 11
 Operaciones en coma flotante: 1378

x = 0.25
x = 0.5
x = 0.75
y = 0.25 y = 0. 5 y = 0.75
6.2447 12.4947 18.7473
12.4947 24.9947 37.4973
18.7473 37.4973 56.2487
Método de Sobrerrelajación.
Factor de relajación: w = 1.2
 Iteraciones: 8
 Operaciones en coma flotante: 1802

x = 0.25
x = 0.5
x = 0.75
y = 0.25 y = 0. 5 y = 0.75
6.2514 12.5008 18.7502
12.5008 25.0003 37.5002
18.7502 37.5002 56.2500
Algoritmos iterativos por bloques

Iteración por bloques fila
Para j = 1, 2, … , m-1, resolver el sistema
- a 2 u(i -k1) , j +  2 + 2a 2  u(i ,kj) - a 2 u(i +k1) , j  u(i ,kj)-1 + u(i ,kj-+11) 


i  1, 2, ..., n 

Iteración por bloques columna
 Método implícito de direcciones alternadas

Método de Direcciones Alternadas.
Iteraciones: 5
 Operaciones en coma flotante: 1468

x = 0.25
x = 0.5
x = 0.75
y = 0.25 y = 0. 5 y = 0.75
6.2490 12.4991 18.7497
12.4989 24.9990 37.4996
18.7495 37.4995 56.2498
Errores máximos.

Solución:
u(x,y) = x·y
Método
Jacobi
Gauss-Seidel
Sobrerrelajación w = 1.2
Direcciones alternadas
y = 0.25
6.2500
x = 0.25
12.5000
x = 0.5
x = 0.75 18.7500
Iteraciones
8
12
8
5
y = 0. 5
12.5000
25.0000
37.5000
Operaciones
1142
1378
1802
1468
y = 0.75
18.7500
37.5000
56.2500
Error máximo
0.0061
0.0053
0.0014
0.0011
FIN