Mtodos Numricos para Ecuaciones en Derivadas Parciales Clculo

Métodos Numéricos para Ecuaciones en Derivadas Parciales Cálculo Numérico Práctica 4

Métodos numéricos para Ecuaciones en Derivadas Parciales n Resolución de sistemas lineales tridiagonales u Factorización de una matriz tridiagonal u Resolución de un sistema mediante LU n Ecuación de Ondas u Método implícito n Ecuación de Laplace u Método de sobrerrelajación por filas

Sistemas lineales tridiagonales A = L·U

Algoritmo de Factorización function [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

Resolución de un sistema mediante LU. Archivo croutlu. m n Eliminación y(1)=d(1)/l(1); for i=2: n y(i)=(d(i)-c(i-1)*y(i-1))/l(i); end n Sustitución regresiva x(n)=y(n); for i=n-1: 1 x(i)=y(i)-u(i)*x(i+1); end

Introducción n n EDP de orden 2, lineales de coeficientes constantes. Auxx+Buxy+Cuyy+Dux+Euy+Fu=G u Ecuación de Ondas utt c 2 uxx = 0 u Ecuación del Calor ut cuxx = 0, c>0 u Ecuación de Laplace uxx + uyy = 0 Condiciones iniciales y de contorno

Ecuación de Ondas n Ecuación de Ondas utt = c²uxx , 0 < x < L, t > 0 n Condiciones iniciales u(x, 0) = f(x) ut(x, 0) = g(x) n Condiciones de contorno u(0, t) = l(t) u(L, t) = r(t) n Ecuación en diferencias finitas

Ejemplo n Ecuación: utt = c²uxx , 0 , < x < L, t > 0, c = 1, L=T=4 n Condiciones iniciales u(x, 0) = 2 |x |x 2|, ut(x, 0) = 0 n Condiciones de contorno u(0, t) = 0, u(L, t) = 0 n Discretización nx=4, nt=8 0 L

Método implícito n Paso 0º: Condición inicial 1ª ui, 0 = fi n Paso 1º: Condición inicial segunda ui, 1 = a = 2 (fi 1+fi+1)/2 + (1 a )/2 + (1 2)fi + kgi n Pasos siguientes: ecuación en diferencias (1+a (1+ 2)ui, j+1 a 2(ui+1, j+1 + ui 1, j+1)/2 = 2 ui, j + a + 2(ui+1, j 1 + ui 1, j 1)/2 (1+a (1+ 2)ui, j 1

Paso del método implícito n Truco ecuación implícita a 2( ui 1, j 1 + ui 1, j+1)/4 + (1 + a + (1 + 2)(ui, j 1 + ui, j+1)/2 a 2(ui+1, j 1 + ui+1, j+1)/4 = ui, j . n Sistema tridiagonal Aw = v, v = (u 1, j, u 2, j, . . . , unx 1, j)' ui, j+1 = wi ui, j 1 n Factorización LU Lz = v Uw = z

Algoritmo: parámetros de entrada nx, h: n nt, k: n c: n f, g: n n hl, hr: nº de intervalos y paso espacial nº de intervalos y paso temporal coeficiente de la ecuación vectores columna (nx+1, 1) de las condiciones iniciales en los nodos con t = 0 vectores fila (1, nt+1)de las condiciones de contorno en los nodos con x = 0 y x = L, resp.

Ejemplo: parámetros de entrada n nx = 4; h=1; x = 0: h: 4; n nt = 8; k =. 5; n c = 1; n f = 2 - abs(2 -x); n g = zeros(size(x)); n hl = zeros(1, nt+1); n hr = zeros(1, nt+1); x = x(: );

Algoritmo: Preparación n a 2 = (c*k/h)^2; % Parámetro de Courant n U = zeros(nx+1, nt+1); n Condiciones de contorno % Solución U(1, : ) = hl; U(nx+1, : ) = hr; n Rangos auxiliares de índices L = 1: nx-1; C = 2: nx; R = 3: nx+1;

Algoritmo: pasos iniciales n Condición inicial sobre la función (paso 0) U(: , 1) = f; n Condición inicial sobre la derivada (paso 1) U(C, 2) = a 2*(f(L) + f(R))/2 +. . . (1 -a 2)*f(C) + k*g(C);

Algoritmo: construcción y factorización de la matriz n Diagonal principal dp = (1+a 2)/2*ones(1, nx-1); n Subdiagonal y superdiagonal ds = -a 2/4*ones(1, nx-2); n Factorización LU [c, l, u]=clu(ds, dp, ds);

Algoritmo: paso general for j = 2: nt % Términos independientes b(1) = a 2/4*(hl(j-1)+hl(j+1)); b(nx-1) = a 2/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

Ejemplo

Ecuaciones elípticas n n Ecuación de Laplace uxx + uyy = 0, 0 < x < a, 0 < y <b Condiciones de contorno u(x, 0) = f 0(x), u(x, b) = f 1(x) u(0, y) = g 0(y), u(a, y) = g 1(y) n Discretización

Ecuación de Laplace n Ecuación en diferencias: a=k/h a 2(ui-1, j + ui+1, j) + ui, j-1 + ui, j+1 2(a 2( 2+1)ui, j = 0 n Matriz del sistema: grande , dispersa n Caso h = k : ui-1, j + ui+1, j + ui, j-1 + ui, j+1 = 4 ui, j

Algoritmos iterativos por bloques Iteración por bloques columna Para j = 1, 2, … , ny 1, resolver el sistema n Iteración por bloques fila n Método implícito de direcciones alternadas n

Ecuación de Laplace. Ejemplo Ecuación uxx+ uyy=0, 0 < x < 1, 0 < y < 1 n Condiciones de contorno u(x, 0) = 0, u (x, 1) = 100 x u(0, y) = 0, u(1, y) = 100 y n Discretización h = 0. 125, k = 0. 25 n

Algoritmo: parámetros de entrada alfa: n f 0, f 1: paso y / paso x vectores columna (nx+1, 1) de las condiciones de contorno en los nodos con y = 0 e y = b, resp. n g 0, g 1: vectores fila (1, ny+1) de las condiciones de contorno en los nodos con x = 0 y x = a, resp. n tol: condición de convergencia n maxiter: tope de iteraciones. n

Ejemplo: parámetros de entrada n h=. 125; x = 0: h: 1; x = x(: ); n k=. 25; n alfa = k/h; n f 0 = zeros(size(x)); n f 1 = 100*x; n g 0 = zeros(size(y)); n g 1 = 100*y; n tol = 5 e-2; n maxiter = 50; y = 0: k: 1;

Algoritmo: Preparación n a 2 = alfa^2; b 2 = 2*(1+a 2); n m = length(f 0); n = length(g 0); n Estimación inicial U = zeros(n, m); n Condiciones de contorno U(: , 1) = f 0; U(: , n) = f 1; U(1, : ) = g 0; U(m, : ) = g 1; 1 g 0 n 1 f 0 f 1 m g 1

Algoritmo: Construcción y factorización de la matriz n Diagonal principal dp = b 2*ones(1, m-2); n Subdiagonal y superdiagonal ds = -a 2*ones(1, m-3); n Factorización LU [c, l, u]=clu(ds, dp, ds);

Algoritmo: relajación por columnas for j = 2: n-1 % Términos independientes b = U(2: m-1, j-1) + U(2: m-1, j+1); b(1) = b(1) + a 2*g 0(j); b(m-2) = b(m-2) + a 2* g 1(j); % Resolucion del sistema U(2: m-1, j) = croutlu(c, l, u, b)'; end

Algoritmo: iteraciones incr = tol + 1; iter = 0; while incr > tol & iter < maxiter Actualizar U por columnas Calcular incr Incrementar iter end

Ejemplo

F I N