Capítulo 3: Métodos de Diferencias Finitas para PDEs Elípticas en 2D
3.1 Introducción
En este capítulo, discutimos métodos de diferencias finitas (FDM) para resolver ecuaciones
diferenciales parciales (PDE) elípticas en dos dimensiones. Estas PDE son fundamentales en
muchas aplicaciones físicas y de ingeniería, incluyendo la electrostática, la mecánica de
fluidos y la transferencia de calor.
3.2 Formulación del Problema
Considere la PDE elíptica general en 2D:
∂²u/∂x² + ∂²u/∂y² = f(x, y)
con las condiciones de contorno adecuadas en el dominio Ω.
• CODIGO MATLAB:
% Parámetros
Lx = 2; % Longitud en x
Ly = 2; % Longitud en y
Nx = 50; % Número de puntos en x
Ny = 50; % Número de puntos en y
% Crear la malla
x = linspace(0, Lx, Nx);
y = linspace(0, Ly, Ny);
dx = x(2) - x(1);
dy = y(2) - y(1);
% Función f(x, y)
f = @(x, y) sin(pi*x).*sin(pi*y);
% Matriz del sistema
A = sparse(Nx*Ny, Nx*Ny);
b = zeros(Nx*Ny, 1);
% Llenar la matriz A y el vector b
for i = 2:Nx-1
for j = 2:Ny-1
k = (j-1)*Nx + i; % índice del punto (i, j)
A(k, k) = -2/dx^2 - 2/dy^2;
A(k, k-1) = 1/dx^2;
A(k, k+1) = 1/dx^2;
A(k, k-Nx) = 1/dy^2;
A(k, k+Nx) = 1/dy^2;
b(k) = f(x(i), y(j));
end
end
% Condiciones de contorno de Dirichlet
for i = 1:Nx
for j = 1:Ny
k = (j-1)*Nx + i; % índice del punto (i, j)
if i == 1 || i == Nx || j == 1 || j == Ny
A(k, :) = 0;
A(k, k) = 1;
b(k) = 0; % Puedes cambiar esto a otra función si tienes una
condición diferente
end
end
end
% Resolver el sistema
u = A\b;
% Convertir el vector solución a una matriz
U = reshape(u, Nx, Ny);
% Graficar la solución
[X, Y] = meshgrid(x, y);
figure;
surf(X, Y, U');
xlabel('x');
ylabel('y');
zlabel('u');
title('Solución de la PDE elíptica');
• RESULTADO:
3.3 Discretización Espacial
Para aplicar el método de diferencias finitas, dividimos el dominio en una cuadrícula de
puntos con espaciado h en ambas direcciones x y y. La derivada segunda en x se puede
aproximar usando la diferencia central:
∂²u/∂x² ≈ (u_{i+1,j} - 2u_{i,j} + u_{i-1,j})/h²
y de manera similar para la derivada en y.
o CÓDIGO MATLAB (SCRIPT):
% Script para resolver una PDE elíptica en 2D usando diferencias finitas
% Parámetros
Lx = 2; % Longitud en x
Ly = 2; % Longitud en y
Nx = 50; % Número de puntos en x
Ny = 50; % Número de puntos en y
% Crear la malla
x = linspace(0, Lx, Nx);
y = linspace(0, Ly, Ny);
h = x(2) - x(1); % Asumimos que dx = dy = h
% Función f(x, y)
f = @(x, y) sin(pi*x).*sin(pi*y);
% Matriz del sistema
A = sparse(Nx*Ny, Nx*Ny);
b = zeros(Nx*Ny, 1);
% Llenar la matriz A y el vector b usando diferencias centrales
for i = 2:Nx-1
for j = 2:Ny-1
k = (j-1)*Nx + i; % índice del punto (i, j)
A(k, k) = -4/h^2;
A(k, k-1) = 1/h^2;
A(k, k+1) = 1/h^2;
A(k, k-Nx) = 1/h^2;
A(k, k+Nx) = 1/h^2;
b(k) = f(x(i), y(j));
end
end
% Condiciones de contorno de Dirichlet
for i = 1:Nx
for j = 1:Ny
k = (j-1)*Nx + i; % índice del punto (i, j)
if i == 1 || i == Nx || j == 1 || j == Ny
A(k, :) = 0;
A(k, k) = 1;
b(k) = 0; % Puedes cambiar esto a otra función si tienes una
condición diferente
end
end
end
% Resolver el sistema
u = A\b;
% Convertir el vector solución a una matriz
U = reshape(u, Nx, Ny);
% Graficar la solución
[X, Y] = meshgrid(x, y);
figure;
surf(X, Y, U');
xlabel('x');
ylabel('y');
zlabel('u');
title('Solución de la PDE elíptica usando diferencias finitas');
• CÓDIGO DE MATLAB (TIPO FUNCIÓN):
function U = solvePDE_discretizacion_function(Lx, Ly, Nx, Ny, f)
% Función para resolver una PDE elíptica en 2D usando diferencias
finitas
% Crear la malla
x = linspace(0, Lx, Nx);
y = linspace(0, Ly, Ny);
h = x(2) - x(1); % Asumimos que dx = dy = h
% Matriz del sistema
A = sparse(Nx*Ny, Nx*Ny);
b = zeros(Nx*Ny, 1);
% Llenar la matriz A y el vector b usando diferencias centrales
for i = 2:Nx-1
for j = 2:Ny-1
k = (j-1)*Nx + i; % índice del punto (i, j)
A(k, k) = -4/h^2;
A(k, k-1) = 1/h^2;
A(k, k+1) = 1/h^2;
A(k, k-Nx) = 1/h^2;
A(k, k+Nx) = 1/h^2;
b(k) = feval(f, x(i), y(j));
end
end
% Condiciones de contorno de Dirichlet
for i = 1:Nx
for j = 1:Ny
k = (j-1)*Nx + i; % índice del punto (i, j)
if i == 1 || i == Nx || j == 1 || j == Ny
A(k, :) = 0;
A(k, k) = 1;
b(k) = 0; % Puedes cambiar esto a otra función si tienes
una condición diferente
end
end
end
% Resolver el sistema
u = A\b;
% Convertir el vector solución a una matriz
U = reshape(u, Nx, Ny);
% Graficar la solución
[X, Y] = meshgrid(x, y);
figure;
surf(X, Y, U');
xlabel('x');
ylabel('y');
zlabel('u');
title('Solución de la PDE elíptica usando diferencias finitas');
end
• RESULTADOS:
3.4 Sistema de Ecuaciones Lineales
La discretización de la PDE elíptica en todos los puntos de la cuadrícula produce un sistema
de ecuaciones lineales de la forma A u = b, donde A es una matriz dispersa que depende de
la discretización y las condiciones de contorno, u es el vector de valores de la solución en los
puntos de la cuadrícula, y b es el vector que contiene los valores de la función f y las
condiciones de contorno.
CÓDIGO MATLAB (TIPO FUNCION)
function U = solvePDE_system_function(Lx, Ly, Nx, Ny, f)
% Función para resolver una PDE elíptica en 2D usando diferencias finitas
% y resolver el sistema de ecuaciones lineales resultante
% Crear la malla
x = linspace(0, Lx, Nx);
y = linspace(0, Ly, Ny);
h = x(2) - x(1); % Asumimos que dx = dy = h
% Matriz del sistema
A = sparse(Nx*Ny, Nx*Ny);
b = zeros(Nx*Ny, 1);
% Llenar la matriz A y el vector b usando diferencias centrales
for i = 2:Nx-1
for j = 2:Ny-1
k = (j-1)*Nx + i; % índice del punto (i, j)
A(k, k) = -4/h^2;
A(k, k-1) = 1/h^2;
A(k, k+1) = 1/h^2;
A(k, k-Nx) = 1/h^2;
A(k, k+Nx) = 1/h^2;
b(k) = feval(f, x(i), y(j));
end
end
% Condiciones de contorno de Dirichlet
for i = 1:Nx
for j = 1:Ny
k = (j-1)*Nx + i; % índice del punto (i, j)
if i == 1 || i == Nx || j == 1 || j == Ny
A(k, :) = 0;
A(k, k) = 1;
b(k) = 0; % Puedes cambiar esto a otra función si tienes una
condición diferente
end
end
end
% Resolver el sistema
u = A\b;
% Convertir el vector solución a una matriz
U = reshape(u, Nx, Ny);
% Graficar la solución
[X, Y] = meshgrid(x, y);
figure;
surf(X, Y, U');
xlabel('x');
ylabel('y');
zlabel('u');
title('Sistema de Ecuaciones Lineales: Solución de la PDE elíptica en
2D');
end
% Parámetros
Lx = 2; % Longitud en x
Ly = 2; % Longitud en y
Nx = 50; % Número de puntos en x
Ny = 50; % Número de puntos en y
% Función f(x, y)
f = @(x, y) sin(pi*x).*sin(pi*y);
% Llamar a la función
U = solvePDE_system_function(Lx, Ly, Nx, Ny, f);
• RESULKTADO:
3.5 Métodos Iterativos
Para resolver el sistema de ecuaciones lineales, a menudo se usan métodos iterativos como
el método de Jacobi, Gauss-Seidel y SOR (Successive Over-Relaxation). Estos métodos son
preferibles a los métodos directos para sistemas grandes debido a su menor requerimiento
de memoria y capacidad para manejar matrices dispersas.
CODIGO DE MATLAB MÉTODO JACOBI:
% Script para resolver una PDE elíptica en 2D usando el método de Jacobi
% Parámetros
Lx = 2; % Longitud en x
Ly = 2; % Longitud en y
Nx = 50; % Número de puntos en x
Ny = 50; % Número de puntos en y
maxIter = 1000; % Número máximo de iteraciones
tol = 1e-6; % Tolerancia para la convergencia
% Crear la malla
x = linspace(0, Lx, Nx);
y = linspace(0, Ly, Ny);
h = x(2) - x(1); % Asumimos que dx = dy = h
% Función f(x, y)
f = @(x, y) sin(pi*x).*sin(pi*y);
% Inicializar la matriz del sistema y el vector b
A = sparse(Nx*Ny, Nx*Ny);
b = zeros(Nx*Ny, 1);
% Llenar la matriz A y el vector b usando diferencias centrales
for i = 2:Nx-1
for j = 2:Ny-1
k = (j-1)*Nx + i; % índice del punto (i, j)
A(k, k) = -4/h^2;
A(k, k-1) = 1/h^2;
A(k, k+1) = 1/h^2;
A(k, k-Nx) = 1/h^2;
A(k, k+Nx) = 1/h^2;
b(k) = f(x(i), y(j));
end
end
% Condiciones de contorno de Dirichlet
for i = 1:Nx
for j = 1:Ny
k = (j-1)*Nx + i; % índice del punto (i, j)
if i == 1 || i == Nx || j == 1 || j == Ny
A(k, :) = 0;
A(k, k) = 1;
b(k) = 0; % Puedes cambiar esto a otra función si tienes una
condición diferente
end
end
end
% Inicializar la solución
u = zeros(Nx*Ny, 1);
u_new = u;
% Método de Jacobi
for iter = 1:maxIter
for i = 2:Nx-1
for j = 2:Ny-1
k = (j-1)*Nx + i; % índice del punto (i, j)
u_new(k) = (b(k) - A(k, k-1)*u(k-1) - A(k, k+1)*u(k+1) - A(k, k-
Nx)*u(k-Nx) - A(k, k+Nx)*u(k+Nx)) / A(k, k);
end
end
% Verificar convergencia
if norm(u_new - u, inf) < tol
break;
end
u = u_new;
end
% Convertir el vector solución a una matriz
U = reshape(u, Nx, Ny);
% Graficar la solución
[X, Y] = meshgrid(x, y);
figure;
surf(X, Y, U');
xlabel('x');
ylabel('y');
zlabel('u');
title('Método de Jacobi: Solución de la PDE elíptica en 2D');
• RESULTADO:
CÓDIGO MATLAB MÉTODO Método de Gauss-Seidel:
% Script para resolver una PDE elíptica en 2D usando el método de Gauss-Seidel
% Parámetros
Lx = 2; % Longitud en x
Ly = 2; % Longitud en y
Nx = 50; % Número de puntos en x
Ny = 50; % Número de puntos en y
maxIter = 1000; % Número máximo de iteraciones
tol = 1e-6; % Tolerancia para la convergencia
% Crear la malla
x = linspace(0, Lx, Nx);
y = linspace(0, Ly, Ny);
h = x(2) - x(1); % Asumimos que dx = dy = h
% Función f(x, y)
f = @(x, y) sin(pi*x).*sin(pi*y);
% Inicializar la matriz del sistema y el vector b
A = sparse(Nx*Ny, Nx*Ny);
b = zeros(Nx*Ny, 1);
% Llenar la matriz A y el vector b usando diferencias centrales
for i = 2:Nx-1
for j = 2:Ny-1
k = (j-1)*Nx + i; % índice del punto (i, j)
A(k, k) = -4/h^2;
A(k, k-1) = 1/h^2;
A(k, k+1) = 1/h^2;
A(k, k-Nx) = 1/h^2;
A(k, k+Nx) = 1/h^2;
b(k) = f(x(i), y(j));
end
end
% Condiciones de contorno de Dirichlet
for i = 1:Nx
for j = 1:Ny
k = (j-1)*Nx + i; % índice del punto (i, j)
if i == 1 || i == Nx || j == 1 || j == Ny
A(k, :) = 0;
A(k, k) = 1;
b(k) = 0; % Puedes cambiar esto a otra función si tienes una
condición diferente
end
end
end
% Inicializar la solución
u = zeros(Nx*Ny, 1);
% Método de Gauss-Seidel
for iter = 1:maxIter
for i = 2:Nx-1
for j = 2:Ny-1
k = (j-1)*Nx + i; % índice del punto (i, j)
u(k) = (b(k) - A(k, k-1)*u(k-1) - A(k, k+1)*u(k+1) - A(k, k-
Nx)*u(k-Nx) - A(k, k+Nx)*u(k+Nx)) / A(k, k);
end
end
% Verificar convergencia
if norm(A*u - b, inf) < tol
break;
end
end
% Convertir el vector solución a una matriz
U = reshape(u, Nx, Ny);
% Graficar la solución
[X, Y] = meshgrid(x, y);
figure;
surf(X, Y, U');
xlabel('x');
ylabel('y');
zlabel('u');
title('Método de Gauss-Seidel: Solución de la PDE elíptica en 2D');
• RESULTADO:
CÓDIGO MATLAB MÉTODO Método SOR:
% Script para resolver una PDE elíptica en 2D usando el método de SOR
% Parámetros
Lx = 2; % Longitud en x
Ly = 2; % Longitud en y
Nx = 50; % Número de puntos en x
Ny = 50; % Número de puntos en y
maxIter = 1000; % Número máximo de iteraciones
tol = 1e-6; % Tolerancia para la convergencia
omega = 1.5; % Factor de relajación
% Crear la malla
x = linspace(0, Lx, Nx);
y = linspace(0, Ly, Ny);
h = x(2) - x(1); % Asumimos que dx = dy = h
% Función f(x, y)
f = @(x, y) sin(pi*x).*sin(pi*y);
% Inicializar la matriz del sistema y el vector b
A = sparse(Nx*Ny, Nx*Ny);
b = zeros(Nx*Ny, 1);
% Llenar la matriz A y el vector b usando diferencias centrales
for i = 2:Nx-1
for j = 2:Ny-1
k = (j-1)*Nx + i; % índice del punto (i, j)
A(k, k) = -4/h^2;
A(k, k-1) = 1/h^2;
A(k, k+1) = 1/h^2;
A(k, k-Nx) = 1/h^2;
A(k, k+Nx) = 1/h^2;
b(k) = f(x(i), y(j));
end
end
% Condiciones de contorno de Dirichlet
for i = 1:Nx
for j = 1:Ny
k = (j-1)*Nx + i; % índice del punto (i, j)
if i == 1 || i == Nx || j == 1 || j == Ny
A(k, :) = 0;
A(k, k) = 1;
b(k) = 0; % Puedes cambiar esto a otra función si tienes una
condición diferente
end
end
end
% Inicializar la solución
u = zeros(Nx*Ny, 1);
% Método SOR
for iter = 1:maxIter
u_old = u; % Guardar la solución anterior
for i = 2:Nx-1
for j = 2:Ny-1
k = (j-1)*Nx + i; % índice del punto (i, j)
u(k) = (1-omega)*u(k) + omega * (b(k) - A(k, k-1)*u(k-1) - A(k,
k+1)*u(k+1) - A(k, k-Nx)*u(k-Nx) - A(k, k+Nx)*u(k+Nx)) / A(k, k);
end
end
% Verificar convergencia
if norm(u - u_old, inf) < tol
break;
end
end
% Convertir el vector solución a una matriz
U = reshape(u, Nx, Ny);
% Graficar la solución
[X, Y] = meshgrid(x, y);
figure;
surf(X, Y, U');
xlabel('x');
ylabel('y');
zlabel('u');
title('Método SOR: Solución de la PDE elíptica en 2D');
• RESULTADO:
3.6 Análisis de Error
El análisis de error de los métodos de diferencias finitas incluye la consistencia, estabilidad
y convergencia del método. La consistencia se refiere a cómo de bien la discretización de
diferencias finitas aproxima la PDE original. La estabilidad asegura que los errores no
crecen exponencialmente en el proceso iterativo. La convergencia asegura que la solución
de diferencias finitas converge a la solución exacta de la PDE cuando el espaciado de la
cuadrícula h tiende a cero.
CÓDIGO MATLAB PARA ANÁLISIS DE CONSISTENCIA:
% Script para analizar la consistencia del método de diferencias finitas
% Parámetros
Lx = 2; % Longitud en x
Ly = 2; % Longitud en y
Nx = 50; % Número de puntos en x
Ny = 50; % Número de puntos en y
h = Lx / (Nx - 1); % Espaciado en x y y
x = linspace(0, Lx, Nx);
y = linspace(0, Ly, Ny);
% Función exacta (solución analítica conocida)
u_exact = @(x, y) sin(pi*x).*sin(pi*y);
% Crear la malla
[X, Y] = meshgrid(x, y);
% Solución exacta
U_exact = u_exact(X, Y);
% Resolver PDE usando diferencias finitas (ejemplo usando método de Jacobi)
% Inicialización
u = zeros(Nx, Ny);
u_new = u;
maxIter = 1000;
tol = 1e-6;
for iter = 1:maxIter
for i = 2:Nx-1
for j = 2:Ny-1
u_new(i, j) = (u(i+1, j) + u(i-1, j) + u(i, j+1) + u(i, j-1)) / 4;
end
end
if norm(u_new - u, inf) < tol
break;
end
u = u_new;
end
% Errores
error = abs(U_exact - u);
max_error = max(max(error));
% Graficar resultados
figure;
subplot(1, 2, 1);
surf(X, Y, U_exact');
title('Solución Exacta');
xlabel('x');
ylabel('y');
zlabel('u_{exact}');
subplot(1, 2, 2);
surf(X, Y, u');
title(['Solución Numérica, Error Máximo = ', num2str(max_error)]);
xlabel('x');
ylabel('y');
zlabel('u_{num}');
• RESULTADO:
CODIGO MATLAB ANÁLISIS DE ESTABILIDAD:
% Script para analizar la estabilidad del método de diferencias finitas
% Parámetros
Lx = 2; % Longitud en x
Ly = 2; % Longitud en y
Nx = 50; % Número de puntos en x
Ny = 50; % Número de puntos en y
h = Lx / (Nx - 1); % Espaciado en x y y
x = linspace(0, Lx, Nx);
y = linspace(0, Ly, Ny);
% Crear la malla
[X, Y] = meshgrid(x, y);
% Función f(x, y)
f = @(x, y) sin(pi*x).*sin(pi*y);
% Inicializar solución
u = zeros(Nx, Ny);
u_new = u;
% Parámetros del método
maxIter = 1000;
tol = 1e-6;
% Método de diferencias finitas (ejemplo usando Jacobi)
error_norms = zeros(maxIter, 1);
for iter = 1:maxIter
% Actualizar la solución
for i = 2:Nx-1
for j = 2:Ny-1
u_new(i, j) = (u(i+1, j) + u(i-1, j) + u(i, j+1) + u(i, j-1)) / 4;
end
end
% Calcular el error
error = u_new - u;
error_norms(iter) = norm(error, 'inf');
% Verificar convergencia
if error_norms(iter) < tol
error_norms = error_norms(1:iter);
break;
end
% Actualizar la solución
u = u_new;
end
% Graficar el análisis de estabilidad
figure;
plot(error_norms, '-o');
xlabel('Iteración');
ylabel('Norma del Error');
title('Análisis de Estabilidad: Norma del Error en Función de Iteraciones');
• RESULTADOS:
CÓDIGO MATLAB PARA ANÁLISIS DE CONVERGENCIA:
% Script para analizar la convergencia del método de diferencias finitas
% Parámetros
Lx = 2; % Longitud en x
Ly = 2; % Longitud en y
h_values = [1/10, 1/20, 1/40, 1/80]; % Espaciado de malla (h)
errors = zeros(length(h_values), 1);
% Función exacta (solución analítica conocida)
u_exact = @(x, y) sin(pi*x).*sin(pi*y);
for idx = 1:length(h_values)
h = h_values(idx);
Nx = round(Lx / h) + 1;
Ny = round(Ly / h) + 1;
x = linspace(0, Lx, Nx);
y = linspace(0, Ly, Ny);
[X, Y] = meshgrid(x, y);
% Inicializar solución
u = zeros(Nx, Ny);
u_new = u;
% Método de diferencias finitas (ejemplo usando Jacobi)
maxIter = 1000;
tol = 1e-6;
for iter = 1:maxIter
for i = 2:Nx-1
for j = 2:Ny-1
u_new(i, j) = (u(i+1, j) + u(i-1, j) + u(i, j+1) + u(i, j-1))
/ 4;
end
end
if norm(u_new - u, inf) < tol
break;
end
u = u_new;
end
% Error
U_exact = u_exact(X, Y);
error = abs(U_exact - u);
errors(idx) = max(max(error));
end
% Graficar el análisis de convergencia
figure;
loglog(h_values, errors, '-o');
xlabel('Tamaño de Malla (h)');
ylabel('Error Máximo');
title('Análisis de Convergencia: Error Máximo en Función del Tamaño de
Malla');
• RESULTADOS:
3.7 Ejemplos Numéricos
Se presentan varios ejemplos numéricos para ilustrar la aplicación de los métodos de
diferencias finitas a problemas de valores de contorno en 2D. Estos ejemplos incluyen la
solución de la ecuación de Poisson y la ecuación de Laplace con diferentes condiciones de
contorno.
CÓDIGO MATLAB PARA LA SOLUCIÓN DE LA ECUACIÓN DE POISSON:
% Script para resolver la ecuación de Poisson en 2D usando diferencias finitas
% Parámetros
Lx = 1; % Longitud en x
Ly = 1; % Longitud en y
Nx = 50; % Número de puntos en x
Ny = 50; % Número de puntos en y
h = Lx / (Nx - 1); % Espaciado en x y y
x = linspace(0, Lx, Nx);
y = linspace(0, Ly, Ny);
% Crear la malla
[X, Y] = meshgrid(x, y);
% Función f(x, y) para la ecuación de Poisson
f = @(x, y) -2*pi^2 * sin(pi*x) * sin(pi*y);
% Inicializar la matriz del sistema y el vector b
A = sparse(Nx*Ny, Nx*Ny);
b = zeros(Nx*Ny, 1);
% Llenar la matriz A y el vector b usando diferencias centrales
for i = 2:Nx-1
for j = 2:Ny-1
k = (j-1)*Nx + i; % índice del punto (i, j)
A(k, k) = -4/h^2;
A(k, k-1) = 1/h^2;
A(k, k+1) = 1/h^2;
A(k, k-Nx) = 1/h^2;
A(k, k+Nx) = 1/h^2;
b(k) = f(x(i), y(j));
end
end
% Condiciones de contorno de Dirichlet (u = 0 en el borde)
for i = 1:Nx
for j = 1:Ny
k = (j-1)*Nx + i; % índice del punto (i, j)
if i == 1 || i == Nx || j == 1 || j == Ny
A(k, :) = 0;
A(k, k) = 1;
b(k) = 0;
end
end
end
% Resolver el sistema lineal
u = A \ b;
% Convertir el vector solución a una matriz
U = reshape(u, Nx, Ny);
% Graficar la solución
figure;
surf(X, Y, U');
xlabel('x');
ylabel('y');
zlabel('u');
title('Solución de la Ecuación de Poisson en 2D');
• RESULTADO:
CÓDIGO MATLAB PARA LA SOLUCIÓN DE LA ECUACIÓN DE LAPLACE CON CONDICIONES
DE CONTORNO DE DIRICHLET:
% Script para resolver la ecuación de Laplace en 2D usando diferencias finitas
% Parámetros
Lx = 1; % Longitud en x
Ly = 1; % Longitud en y
Nx = 50; % Número de puntos en x
Ny = 50; % Número de puntos en y
h = Lx / (Nx - 1); % Espaciado en x y y
x = linspace(0, Lx, Nx);
y = linspace(0, Ly, Ny);
% Crear la malla
[X, Y] = meshgrid(x, y);
% Función para el lado derecho (en la ecuación de Laplace, esto es cero)
f = @(x, y) 0;
% Inicializar la matriz del sistema y el vector b
A = sparse(Nx*Ny, Nx*Ny);
b = zeros(Nx*Ny, 1);
% Llenar la matriz A y el vector b usando diferencias centrales
for i = 2:Nx-1
for j = 2:Ny-1
k = (j-1)*Nx + i; % índice del punto (i, j)
A(k, k) = -4/h^2;
A(k, k-1) = 1/h^2;
A(k, k+1) = 1/h^2;
A(k, k-Nx) = 1/h^2;
A(k, k+Nx) = 1/h^2;
b(k) = f(x(i), y(j));
end
end
% Condiciones de contorno de Dirichlet (por ejemplo, u = x*(1-x)*y*(1-y))
for i = 1:Nx
for j = 1:Ny
k = (j-1)*Nx + i; % índice del punto (i, j)
if i == 1
A(k, :) = 0;
A(k, k) = 1;
b(k) = x(i)*(1-x(i))*y(j)*(1-y(j)); % Contorno en x = 0
elseif i == Nx
A(k, :) = 0;
A(k, k) = 1;
b(k) = x(i)*(1-x(i))*y(j)*(1-y(j)); % Contorno en x = Lx
elseif j == 1
A(k, :) = 0;
A(k, k) = 1;
b(k) = x(i)*(1-x(i))*y(j)*(1-y(j)); % Contorno en y = 0
elseif j == Ny
A(k, :) = 0;
A(k, k) = 1;
b(k) = x(i)*(1-x(i))*y(j)*(1-y(j)); % Contorno en y = Ly
end
end
end
% Resolver el sistema lineal
u = A \ b;
% Convertir el vector solución a una matriz
U = reshape(u, Nx, Ny);
% Graficar la solución
figure;
surf(X, Y, U');
xlabel('x');
ylabel('y');
zlabel('u');
title('Solución de la Ecuación de Laplace en 2D con Condiciones de Contorno de
Dirichlet');
• RESULTADOS:
CÓDIGO DE MATLAB PARA LA SOLUCIÓN DE LA ECUACIÓN DE LA PLACE CON
CONDICIONES DE CONTORNO DE NUEMANN:
% Script para resolver la ecuación de Laplace en 2D usando diferencias finitas
con condiciones de Neumann
% Parámetros
Lx = 1; % Longitud en x
Ly = 1; % Longitud en y
Nx = 50; % Número de puntos en x
Ny = 50; % Número de puntos en y
h = Lx / (Nx - 1); % Espaciado en x y y
x = linspace(0, Lx, Nx);
y = linspace(0, Ly, Ny);
% Crear la malla
[X, Y] = meshgrid(x, y);
% Función para el lado derecho (en la ecuación de Laplace, esto es cero)
f = @(x, y) 0;
% Inicializar la matriz del sistema y el vector b
A = sparse(Nx*Ny, Nx*Ny);
b = zeros(Nx*Ny, 1);
% Llenar la matriz A y el vector b usando diferencias centrales
for i = 2:Nx-1
for j = 2:Ny-1
k = (j-1)*Nx + i; % índice del punto (i, j)
A(k, k) = -4/h^2;
A(k, k-1) = 1/h^2;
A(k, k+1) = 1/h^2;
A(k, k-Nx) = 1/h^2;
A(k, k+Nx) = 1/h^2;
b(k) = f(x(i), y(j));
end
end
% Condiciones de contorno de Neumann (∂u/∂n = 0 en el borde)
for i = 1:Nx
for j = 1:Ny
k = (j-1)*Nx + i; % índice del punto (i, j)
if i == 1 % Condición en x = 0
A(k, k) = 1;
b(k) = 0; % ∂u/∂x = 0
elseif i == Nx % Condición en x = Lx
A(k, k) = 1;
b(k) = 0; % ∂u/∂x = 0
elseif j == 1 % Condición en y = 0
A(k, k) = 1;
b(k) = 0; % ∂u/∂y = 0
elseif j == Ny % Condición en y = Ly
A(k, k) = 1;
b(k) = 0; % ∂u/∂y = 0
end
end
end
% Resolver el sistema lineal
u = A \ b;
% Convertir el vector solución a una matriz
U = reshape(u, Nx, Ny);
% Graficar la solución
figure;
surf(X, Y, U');
xlabel('x');
ylabel('y');
zlabel('u');
title('Solución de la Ecuación de Laplace en 2D con Condiciones de Contorno de
Neumann');
• RESULTADOS:
3.8 Programación de Métodos de Diferencias Finitas en 2D
La implementación computacional de los métodos de diferencias finitas se realiza mediante
la creación de programas que ensamblan la matriz A y el vector b, y luego resuelven el
sistema lineal resultante usando métodos iterativos. Se proporciona un esquema general de
un programa en MATLAB o Python, incluyendo la configuración de la cuadrícula, la
aplicación de las condiciones de contorno y la iteración para resolver el sistema.
CÓDIGO MATLAB Programación de Métodos de Diferencias Finitas en 2D con Método
de Gauss-Seidel:
% Parámetros del dominio y discretización
Lx = 1; % Longitud en x
Ly = 1; % Longitud en y
Nx = 50; % Número de puntos en x
Ny = 50; % Número de puntos en y
h = Lx / (Nx - 1); % Espaciado en x y y
x = linspace(0, Lx, Nx);
y = linspace(0, Ly, Ny);
% Función f(x, y) para la ecuación de Poisson
f = @(x, y) -2*pi^2 * sin(pi*x) * sin(pi*y);
% Inicialización de la matriz A y el vector b
A = sparse(Nx*Ny, Nx*Ny);
b = zeros(Nx*Ny, 1);
% Ensamblado de la matriz A y el vector b usando diferencias finitas
for i = 2:Nx-1
for j = 2:Ny-1
k = (j-1)*Nx + i; % índice del punto (i, j)
A(k, k) = -4/h^2;
A(k, k-1) = 1/h^2;
A(k, k+1) = 1/h^2;
A(k, k-Nx) = 1/h^2;
A(k, k+Nx) = 1/h^2;
b(k) = f(x(i), y(j));
end
end
% Condiciones de contorno de Dirichlet (u = 0 en el borde)
for i = 1:Nx
for j = 1:Ny
k = (j-1)*Nx + i; % índice del punto (i, j)
if i == 1 || i == Nx || j == 1 || j == Ny
A(k, :) = 0;
A(k, k) = 1;
b(k) = 0;
end
end
end
% Método de Gauss-Seidel para resolver el sistema lineal
u0 = zeros(Nx*Ny, 1); % Vector inicial de aproximación
maxIter = 1000; % Número máximo de iteraciones
tol = 1e-6; % Tolerancia de convergencia
u = gauss_seidel_method(A, b, u0, maxIter, tol);
% Convertir el vector solución a una matriz
U = reshape(u, Nx, Ny);
% Graficar la solución
[X, Y] = meshgrid(x, y);
figure;
surf(X, Y, U');
xlabel('x');
ylabel('y');
zlabel('u');
title('Solución de la Ecuación de Poisson en 2D con Condiciones de Contorno de
Dirichlet (Método de Gauss-Seidel)');
% Función para el método de Gauss-Seidel
function u = gauss_seidel_method(A, b, u0, maxIter, tol)
L = tril(A); % Parte triangular inferior de A
U = triu(A, 1); % Parte triangular superior de A
u = u0;
for iter = 1:maxIter
u_old = u;
for i = 1:length(b)
u(i) = (b(i) - U(i, :) * u) / L(i, i);
end
if norm(u - u_old, inf) < tol
fprintf('Gauss-Seidel convergió en %d iteraciones.\n', iter);
break;
end
end
end
• RESULTADOS:
Ejercicios
1. Demostrar la consistencia del esquema: Demuestre que el siguiente esquema es
consistente con la PDE ∂u/∂t + C ∂u/∂x = F:
u_i^{n+1} = u_i^n - Δt/2h (u_{i+1}^n - u_{i-1}^n) + C Δt ( (u_{i+1}^n - u_{i-1}^n) / 2h ) + a Δt
( (u_{i+1}^n - 2u_i^n + u_{i-1}^n) / h² )
Discuta también la estabilidad, en la medida de lo posible.
• SOLUCIÓN EJERCICIO 1:
% Parámetros
C = 1; % Velocidad de propagación
a = 1; % Coeficiente de difusión
h = 0.1; % Paso espacial
dt = 0.01; % Paso temporal
% Discretización espacial y temporal
x = 0:h:1; % Dominio espacial
t = 0:dt:0.1; % Dominio temporal
nx = length(x); % Número de puntos espaciales
nt = length(t); % Número de puntos temporales
% Condiciones iniciales
u = zeros(nx, nt);
u(:,1) = sin(pi*x); % Ejemplo de condición inicial
% Implementación del esquema numérico
for n = 1:nt-1
for i = 2:nx-1
u(i,n+1) = u(i,n) - (dt/(2*h)) * (u(i+1,n) - u(i-1,n)) ...
+ C * dt * ((u(i+1,n) - u(i-1,n)) / (2*h)) ...
+ a * dt * ((u(i+1,n) - 2*u(i,n) + u(i-1,n)) / (h^2));
end
end
% Gráfica de la solución numérica
figure;
surf(t, x, u');
title('Solución numérica');
xlabel('Tiempo');
ylabel('Espacio');
zlabel('u');
shading interp;
% Análisis de consistencia
syms u_i_n u_i_n1 u_i_m1 u_ip1 u_im1 dt h
LHS = u_i_n1 - u_i_n;
RHS = - (dt/(2*h)) * (u_ip1 - u_im1) + C * dt * ((u_ip1 - u_im1) / (2*h)) ...
+ a * dt * ((u_ip1 - 2*u_i_n + u_im1) / (h^2));
% Expandir en series de Taylor alrededor de (i,n)
LHS_taylor = taylor(u_i_n1 - u_i_n, dt);
RHS_taylor = taylor(RHS, [u_i_n, u_ip1, u_im1], 'Order', 3);
disp('LHS Taylor expansion:');
disp(LHS_taylor);
disp('RHS Taylor expansion:');
disp(RHS_taylor);
% Discusión de estabilidad (condición de CFL)
CFL = C * dt / h;
if CFL <= 1
disp('El esquema es condicionalmente estable bajo la condición CFL <= 1');
else
disp('El esquema no es estable bajo la condición CFL > 1');
end
LHS Taylor expansion:
u_i_n1 - u_i_n
RHS Taylor expansion:
(dt*u_im1)/h^2 + (dt*u_ip1)/h^2 - (2*dt*u_i_n)/h^2
• RESLTADO DE EJERCICIO 1:
2. Implementar y probar esquemas: Implemente y pruebe los esquemas de diferencias
finitas centradas y de Lax-Wendroff para la ecuación de onda unidireccional ∂u/∂t + ∂u/∂x
= 0. Suponga que el dominioes -1 ≤ x ≤ 1 y t_final = 1. Pruebe su código para los siguientes
parámetros:
(a) u(t, -1) = 0, y u(0, x) = { (x + 1)² / 4 si -1 ≤ x ≤ 0; (1 - x)² / 4 si 0 ≤ x ≤ 1 }
(b) u(t, -1) = 0, y u(0, x) = 1.
Haga el análisis de refinamiento de la malla en t_final = 1 para el caso (a) donde la solución
exacta está disponible, tomando h = 0.1, 0.05, 0.025, 0.0125. Para el problema (b), utilice h =
0.025. Grafique la solución en t_final = 1 para ambos casos.
SOLUCIÓN:
o Archivo Principal: (main_script.m)
% Parámetros del problema
x_min = -1;
x_max = 1;
t_final = 1;
CFL = 0.5; % Número de Courant-Friedrichs-Lewy
% Discretización del dominio
h_values = [0.1, 0.05, 0.025, 0.0125]; % Diferentes tamaños de malla
para análisis de refinamiento
% Condiciones iniciales caso (a)
initial_condition_a = @(x) (x + 1).^2 / 4 .* (x >= -1 & x <= 0) + (1 -
x).^2 / 4 .* (x > 0 & x <= 1);
% Condiciones iniciales caso (b)
initial_condition_b = @(x) 1;
% Generar gráficos
figure;
% Caso (a)
for i = 1:length(h_values)
h = h_values(i);
x_points = x_min:h:x_max;
u = finite_difference_scheme(x_points, t_final, h, 'lax-wendroff',
initial_condition_a);
subplot(2, 2, i);
plot(x_points, u(:, end), 'DisplayName', sprintf('h = %.4f', h));
title(sprintf('Caso (a) - Lax-Wendroff - h = %.4f', h));
xlabel('x');
ylabel('u');
legend;
grid on;
end
% Refinamiento de malla - Caso (a)
figure;
errors = zeros(length(h_values), 1);
exact_solution = @(x, t) initial_condition_a(x - t);
for i = 1:length(h_values)
h = h_values(i);
x_points = x_min:h:x_max;
u = finite_difference_scheme(x_points, t_final, h, 'lax-wendroff',
initial_condition_a);
exact = exact_solution(x_points, t_final);
errors(i) = max(abs(u(:, end) - exact'));
subplot(2, 2, i);
plot(x_points, u(:, end), 'b', x_points, exact, 'r--');
title(sprintf('Refinamiento de malla - h = %.4f', h));
xlabel('x');
ylabel('u');
legend('Numérica', 'Exacta');
grid on;
end
figure;
loglog(h_values, errors, '-o');
title('Error de refinamiento de malla - Caso (a)');
xlabel('h');
ylabel('Error máximo');
grid on;
% Caso (b)
h = 0.025;
x_points = x_min:h:x_max;
u = finite_difference_scheme(x_points, t_final, h, 'lax-wendroff',
initial_condition_b);
figure;
plot(x_points, u(:, end));
title('Caso (b) - Lax-Wendroff - h = 0.025');
xlabel('x');
ylabel('u');
grid on;
o ARCHIVO DE FUNCIÓN: (finite_difference_scheme.m)
function u = finite_difference_scheme(x, t_final, h, scheme,
initial_condition)
x_min = min(x);
x_max = max(x);
CFL = 0.5; % Número de Courant-Friedrichs-Lewy
k = CFL * h; % Paso de tiempo determinado por CFL
x_points = x_min:h:x_max;
t_points = 0:k:t_final;
u = zeros(length(x_points), length(t_points));
% Condiciones iniciales
u(:, 1) = initial_condition(x_points);
% Esquemas numéricos
for n = 1:length(t_points)-1
for i = 2:length(x_points)-1
if strcmp(scheme, 'centrado')
u(i, n+1) = u(i, n) - CFL/2 * (u(i+1, n) - u(i-1, n));
elseif strcmp(scheme, 'lax-wendroff')
u(i, n+1) = u(i, n) - CFL/2 * (u(i+1, n) - u(i-1, n)) +
...
CFL^2/2 * (u(i+1, n) - 2*u(i, n) + u(i-1,
n));
end
end
% Condición de frontera
u(1, n+1) = 0;
u(end, n+1) = u(end-1, n+1); % Aproximación a la frontera
derecha
end
end
o RESULTADOS EJERCICIO 2:
3. Resolver la ecuación de Burgers: Utilice los esquemas de diferencias finitas centradas y de
Lax-Wendroff para la ecuación de Burgers:
∂u/∂t + u ∂u/∂x = 0
con las mismas condiciones iniciales y de contorno que en el problema 2.
• SOLUCIÓN PROBLEMA 3:
o Archivo principal (Scrip) :
% Parámetros del problema
x_min = -1;
x_max = 1;
t_final = 1;
CFL = 0.5; % Número de Courant-Friedrichs-Lewy
% Discretización del dominio
h_values = [0.1, 0.05, 0.025, 0.0125]; % Diferentes tamaños de malla
para análisis de refinamiento
% Condiciones iniciales caso (a)
initial_condition_a = @(x) (x + 1).^2 / 4 .* (x >= -1 & x <= 0) + (1 -
x).^2 / 4 .* (x > 0 & x <= 1);
% Condiciones iniciales caso (b)
initial_condition_b = @(x) 1;
% Generar gráficos
figure;
% Caso (a)
for i = 1:length(h_values)
h = h_values(i);
x_points = x_min:h:x_max;
u = burgers_scheme(x_points, t_final, h, 'lax-wendroff',
initial_condition_a, CFL);
subplot(2, 2, i);
plot(x_points, u(:, end), 'DisplayName', sprintf('h = %.4f', h));
title(sprintf('Caso (a) - Lax-Wendroff - h = %.4f', h));
xlabel('x');
ylabel('u');
legend;
grid on;
end
% Refinamiento de malla - Caso (a)
figure;
errors = zeros(length(h_values), 1);
exact_solution = @(x, t) initial_condition_a(x - t); % Aproximación de
la solución exacta para la comparación
for i = 1:length(h_values)
h = h_values(i);
x_points = x_min:h:x_max;
u = burgers_scheme(x_points, t_final, h, 'lax-wendroff',
initial_condition_a, CFL);
exact = exact_solution(x_points, t_final);
errors(i) = max(abs(u(:, end) - exact'));
subplot(2, 2, i);
plot(x_points, u(:, end), 'b', x_points, exact, 'r--');
title(sprintf('Refinamiento de malla - h = %.4f', h));
xlabel('x');
ylabel('u');
legend('Numérica', 'Exacta');
grid on;
end
figure;
loglog(h_values, errors, '-o');
title('Error de refinamiento de malla - Caso (a)');
xlabel('h');
ylabel('Error máximo');
grid on;
% Caso (b)
h = 0.025;
x_points = x_min:h:x_max;
u = burgers_scheme(x_points, t_final, h, 'lax-wendroff',
initial_condition_b, CFL);
figure;
plot(x_points, u(:, end));
title('Caso (b) - Lax-Wendroff - h = 0.025');
xlabel('x');
ylabel('u');
grid on;
o Archivo de Función:
function u = burgers_scheme(x, t_final, h, scheme, initial_condition,
CFL)
x_min = min(x);
x_max = max(x);
k = CFL * h; % Paso de tiempo determinado por CFL
x_points = x_min:h:x_max;
t_points = 0:k:t_final;
u = zeros(length(x_points), length(t_points));
% Condiciones iniciales
u(:, 1) = initial_condition(x_points);
% Esquemas numéricos
for n = 1:length(t_points)-1
for i = 2:length(x_points)-1
if strcmp(scheme, 'centrado')
u(i, n+1) = u(i, n) - CFL/2 * (u(i+1, n)^2 - u(i-1,
n)^2) / (2*h);
elseif strcmp(scheme, 'lax-wendroff')
u(i, n+1) = u(i, n) - CFL/2 * (u(i+1, n)^2 - u(i-1,
n)^2) / (2*h) + ...
CFL^2/2 * (u(i+1, n) - 2*u(i, n) + u(i-1,
n));
end
end
% Condición de frontera
u(1, n+1) = 0;
u(end, n+1) = u(end-1, n+1); % Aproximación a la frontera
derecha
end
end
• RESULTADOS EJERCICIO 3:
4. Resolver numéricamente la ecuación de onda: Resuelva numéricamente la siguiente
ecuación de onda usando un método de diferencias finitas:
∂²u/∂t² = c² ∂²u/∂x² + F(x, t)
para 0 < x < 1, con F(x, t) = { x/4 si 0 ≤ x < 0.5; (1 - x)/4 si 0.5 ≤ x < 1 }.
Pruebe su código y orden de convergencia usando un problema que tenga la solución
exacta.
• SOLUCIÓN EJERCICIO 4:
o Archivo principal (wave_solver.m):
% Parámetros del problema
x_min = 0;
x_max = 1;
t_final = 1;
c = 1; % Velocidad de propagación de la onda
CFL = 0.5; % Número de Courant-Friedrichs-Lewy
% Función F(x, t)
F = @(x, t) (x/4) .* (x >= 0 & x < 0.5) + ((1 - x)/4) .* (x >= 0.5 & x
<= 1);
% Discretización del dominio
h_values = [0.1, 0.05, 0.025, 0.0125]; % Diferentes tamaños de malla
para análisis de refinamiento
% Generar gráficos
figure;
% Refinamiento de malla
errors = zeros(length(h_values), 1);
exact_solution = @(x, t) sin(pi*x)*cos(pi*t); % Solución exacta para
comparar
for i = 1:length(h_values)
h = h_values(i);
x_points = x_min:h:x_max;
u = wave_equation_solver(x_points, t_final, h, c, F, CFL, x_min,
x_max);
exact = exact_solution(x_points, t_final);
errors(i) = max(abs(u(:, end) - exact'));
subplot(2, 2, i);
plot(x_points, u(:, end), 'b', x_points, exact, 'r--');
title(sprintf('Refinamiento de malla - h = %.4f', h));
xlabel('x');
ylabel('u');
legend('Numérica', 'Exacta');
grid on;
end
figure;
loglog(h_values, errors, '-o');
title('Error de refinamiento de malla');
xlabel('h');
ylabel('Error máximo');
grid on;
% Prueba con una condición inicial
h = 0.025;
x_points = x_min:h:x_max;
u = wave_equation_solver(x_points, t_final, h, c, F, CFL, x_min, x_max);
figure;
plot(x_points, u(:, end));
title('Solución numérica - h = 0.025');
xlabel('x');
ylabel('u');
grid on;
o Archivo Función (wave_equation_solver.m):
function u = wave_equation_solver(x_points, t_final, h, c, F, CFL,
x_min, x_max)
k = CFL * h / c; % Paso de tiempo determinado por CFL
t_points = 0:k:t_final;
u = zeros(length(x_points), length(t_points));
u_prev = zeros(length(x_points), 1);
% Condiciones iniciales
% u(x,0) = 0 y ∂u/∂t (x,0) = 0
u(:, 1) = 0; % u(x, 0) = 0
u(:, 2) = 0; % ∂u/∂t (x, 0) = 0
% Esquema numérico
for n = 2:length(t_points)-1
for i = 2:length(x_points)-1
u(i, n+1) = 2*(1 - CFL^2)*u(i, n) - u(i, n-1) +
CFL^2*(u(i+1, n) + u(i-1, n)) + k^2*F(x_points(i), t_points(n));
end
% Condiciones de contorno
u(1, n+1) = 0; % Frontera izquierda
u(end, n+1) = 0; % Frontera derecha
end
end
• RESULTADOS DE EJERCICIO 4:
