1 Métodos Numéricos para Ecuaciones en Derivadas ParcialesCálculo Numérico Práctica 4
2 Métodos numéricos para Ecuaciones en Derivadas ParcialesResolución de sistemas lineales tridiagonales Factorización de una matriz tridiagonal Resolución de un sistema mediante LU Ecuación de Ondas Método implícito Ecuación de Laplace Método de sobrerrelajación por filas
3 Sistemas lineales tridiagonalesA = L·U
4 Algoritmo de Factorizaciónfunction [c,l,u]=clu(c,a,b) n = length(a); l(1) = a(1); for i=2:n u(i-1) = b(i-1)/l(i-1); l(i) = a(i) - c(i-1)*u(i-1); end
5 Resolución de un sistema mediante LU. Archivo croutlu.mEliminación y(1)=d(1)/l(1); for i=2:n y(i)=(d(i)-c(i-1)*y(i-1))/l(i); end Sustitución regresiva x(n)=y(n); for i=n-1:-1:1 x(i)=y(i)-u(i)*x(i+1);
6 Auxx+Buxy+Cuyy+Dux+Euy+Fu=GIntroducció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, c>0 Ecuación de Laplace uxx + uyy = 0 Condiciones iniciales y de contorno
7 Ecuación de Ondas Ecuación de Ondas Condiciones inicialesCondiciones de contorno Ecuación en diferencias finitas utt = c²uxx , 0 < x < L, t > 0 u(x, 0) = f(x) ut(x, 0) = g(x) u(0,t) = l(t) u(L,t) = r(t)
8 Ejemplo Ecuación: utt = c²uxx , 0 < x < L, t > 0, c = 1, L=T=4 Condiciones iniciales u(x, 0) = 2|x-2|, ut(x, 0) = 0 Condiciones de contorno u(0,t) = 0, u(L,t) = 0 Discretización nx=4, nt=8 L
9 ui,1 = a2 (fi-1+fi+1)/2 + (1-a2)fi + kgiMétodo implícito Paso 0º: Condición inicial 1ª ui,0 = fi Paso 1º: Condición inicial segunda ui,1 = a2 (fi-1+fi+1)/2 + (1-a2)fi + kgi Pasos siguientes: ecuación en diferencias (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
10 Paso del método implícitoTruco 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 Aw = v, v = (u1,j,u2,j,...,unx-1,j)' tridiagonal ui,j+1 = wi - ui,j-1 Factorización LU Lz = v Uw = z
11 Algoritmo: parámetros de entradanx, h: nº de intervalos y paso espacial nt, k: nº de intervalos y paso temporal c: coeficiente de la ecuación f, g: vectores columna (nx+1,1) de las condiciones iniciales en los nodos con t = 0 hl, hr: vectores fila (1, nt+1)de las condiciones de contorno en los nodos con x = 0 y x = L, resp.
12 Ejemplo: parámetros de entradanx = 4; h=1; x = 0:h:4; x = x(:); nt = 8; k =.5; c = 1; f = 2 - abs(2-x); g = zeros(size(x)); hl = zeros(1, nt+1); hr = zeros(1, nt+1);
13 Algoritmo: Preparacióna2 = (c*k/h)^2; % Parámetro de Courant U = zeros(nx+1,nt+1); % Solución Condiciones de contorno U(1,:) = hl; U(nx+1,:) = hr; Rangos auxiliares de índices L = 1:nx-1; C = 2:nx; R = 3:nx+1;
14 Algoritmo: pasos inicialesCondición inicial sobre la función (paso 0) U(:,1) = f; Condición inicial sobre la derivada (paso 1) U(C,2) = a2*(f(L) + f(R))/ (1-a2)*f(C) + k*g(C);
15 Algoritmo: construcción y factorización de la matrizDiagonal principal dp = (1+a2)/2*ones(1, nx-1); Subdiagonal y superdiagonal ds = -a2/4*ones(1, nx-2); Factorización LU [c,l,u]=clu(ds,dp,ds);
16 Algoritmo: paso generalfor j = 2:nt % Términos independientes b(1) = a2/4*(hl(j-1)+hl(j+1)); b(nx-1) = a2/4*(hr(j-1)+hr(j+1)); % Resolucion del sistema U(C,j+1) = ... croutlu(c,l,u,U(C,j)+b')'-U(C,j-1); end
17 Ejemplo
18 Ecuaciones elípticas Ecuación de Laplaceuxx + uyy = 0, < x < a, < y
19 Ecuación de Laplace Ecuación en diferencias: a=k/ha2(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
20 Algoritmos iterativos por bloquesIteración por bloques columna Para j = 1, 2, … , ny-1, resolver el sistema Iteración por bloques fila Método implícito de direcciones alternadas
21 Ecuación de Laplace. Ejemplouxx+ uyy=0, 0 < x < 1, < y < 1 Condiciones de contorno u(x, 0) = 0, u (x, 1) = 100x u(0, y) = 0, u(1, y) = 100y Discretización h = 0.125, k = 0.25
22 Algoritmo: parámetros de entradaalfa: paso y / paso x f0, f1: vectores columna (nx+1, 1) de las condiciones de contorno en los nodos con y = 0 e y = b, resp. g0, g1: vectores fila (1, ny+1) de las condiciones de contorno en los nodos con x = 0 y x = a, resp. tol: condición de convergencia maxiter: tope de iteraciones.
23 Ejemplo: parámetros de entradah=.125; x = 0:h:1; x = x(:); k=.25; y = 0:k:1; alfa = k/h; f0 = zeros(size(x)); f1 = 100*x; g0 = zeros(size(y)); g1 = 100*y; tol = 5e-2; maxiter = 50;
24 Algoritmo: Preparacióna2 = alfa^2; b2 = 2*(1+a2); m = length(f0); n = length(g0); Estimación inicial U = zeros(n, m); Condiciones de contorno U(:,1) = f0; U(:,n) = f1; U(1,:) = g0; U(m,:) = g1; g n 1 f0 m f1 g1
25 Algoritmo: Construcción y factorización de la matrizDiagonal principal dp = b2*ones(1, m-2); Subdiagonal y superdiagonal ds = -a2*ones(1, m-3); Factorización LU [c,l,u]=clu(ds,dp,ds);
26 Algoritmo: relajación por columnasfor j = 2:n-1 % Términos independientes b = U(2:m-1, j-1) + U(2:m-1, j+1); b(1) = b(1) + a2*g0(j); b(m-2) = b(m-2) + a2* g1(j); % Resolucion del sistema U(2:m-1, j) = croutlu(c,l,u,b)'; end
27 Algoritmo: iteracionesincr = tol + 1; iter = 0; while incr > tol & iter < maxiter Actualizar U por columnas Calcular incr Incrementar iter end
28 Ejemplo
29 F I N