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
a1
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