El método de diferencias finitas es una de las técnicas numéricas más extendidas para la resolución de ecuaciones en derivadas parciales con condiciones de contorno.

Fundamentos teóricos

Formulación del problema

Se considera la ecuación de Poisson en dos dimensiones:

-∇²u = f(x,y), en el dominio Ω
u = g(x,y), en el contorno ∂Ω

donde ∇² representa el operador laplaciano, f el término fuente y g las condiciones de frontera.

Principio discretización

La idea fundamental consiste en aproximar derivadas mediante cocientes diferenciales, transformando la ecuación diferencial en un sistema algebraico.

Aproximación del operador laplaciano:

∇²u ≈ [u(x+h,y) + u(x-h,y) + u(x,y+h) + u(x,y-h) - 4u(x,y)] / h²

Implementación numérica

Esquema de cinco puntos

Para una malla uniforme, se emplea comúnmente el esquema de cinco puntos:

function [U, X, Y] = resolverEDP2D(f, g, xmin, xmax, ymin, ymax, nx, ny)
    % Resuelve una EDP bidimensional mediante diferencias finitas
    % Entradas:
    %   f: término fuente f(x,y)
    %   g: condición de frontera g(x,y)
    %   xmin, xmax: límites en x
    %   ymin, ymax: límites en y
    %   nx, ny: número de puntos en cada dirección

    dx = (xmax - xmin) / (nx - 1);
    dy = (ymax - ymin) / (ny - 1);

    X = linspace(xmin, xmax, nx);
    Y = linspace(ymin, ymax, ny);
    [Xmalla, Ymalla] = meshgrid(X, Y);

    U = zeros(ny, nx);

    % Establecer condiciones de frontera
    for idx = 1:nx
        U(1, idx) = g(X(idx), Y(1));
        U(end, idx) = g(X(idx), Y(end));
    end
    for jdx = 1:ny
        U(jdx, 1) = g(X(1), Y(jdx));
        U(jdx, end) = g(X(end), Y(jdx));
    end

    % Construir sistema lineal
    Nnodos = (nx-2)*(ny-2);
    A = sparse(Nnodos, Nnodos);
    B = zeros(Nnodos, 1);

    alpha_x = 1/dx^2;
    alpha_y = 1/dy^2;
    beta = 2*(alpha_x + alpha_y);

    for jj = 2:ny-1
        for ii = 2:nx-1
            k = (jj-2)*(nx-2) + (ii-1);

            A(k,k) = beta;

            if ii > 2
                A(k,k-1) = -alpha_x;
            end
            if ii < nx-1
                A(k,k+1) = -alpha_x;
            end
            if jj > 2
                A(k,k-(nx-2)) = -alpha_y;
            end
            if jj < ny-1
                A(k,k+(nx-2)) = -alpha_y;
            end

            B(k) = f(X(ii), Y(jj));

            % Incorporar contribuciones de frontera
            if ii == 2
                B(k) = B(k) + alpha_x * g(X(1), Y(jj));
            end
            if ii == nx-1
                B(k) = B(k) + alpha_x * g(X(end), Y(jj));
            end
            if jj == 2
                B(k) = B(k) + alpha_y * g(X(ii), Y(1));
            end
            if jj == ny-1
                B(k) = B(k) + alpha_y * g(X(ii), Y(end));
            end
        end
    end

    % Resolver sistema
    solucion_interna = A \ B;

    % Asignar solución interior
    for jj = 2:ny-1
        for ii = 2:nx-1
            k = (jj-2)*(nx-2) + (ii-1);
            U(jj,ii) = solucion_interna(k);
        end
    end
end

Método de relajación sucesiva (SOR)

Para problemas de gran escala, los métodos iterativos resultan más eficientes:

function U = resolverSOR(f, g, xmin, xmax, ymin, ymax, nx, ny, omega, maxiter, tol)
    % Resuelve mediante SOR
    dx = (xmax - xmin) / (nx - 1);
    dy = (ymax - ymin) / (ny - 1);
    X = linspace(xmin, xmax, nx);
    Y = linspace(ymin, ymax, ny);

    U = zeros(ny, nx);

    % Condiciones de frontera iniciales
    for idx = 1:nx
        U(1,idx) = g(X(idx), Y(1));
        U(end,idx) = g(X(idx), Y(end));
    end
    for jdx = 1:ny
        U(jdx,1) = g(X(1), Y(jdx));
        U(jdx,end) = g(X(end), Y(jdx));
    end

    alpha_x = 1/dx^2;
    alpha_y = 1/dy^2;
    gamma = 2*(alpha_x + alpha_y);

    for iter = 1:maxiter
        error_max = 0;

        for jj = 2:ny-1
            for ii = 2:nx-1
                residuo = (alpha_x*(U(jj,ii-1) + U(jj,ii+1)) + ...
                           alpha_y*(U(jj-1,ii) + U(jj+1,ii)) - ...
                           f(X(ii), Y(jj))) / gamma - U(jj,ii);
                
                U(jj,ii) = U(jj,ii) + omega * residuo;
                error_max = max(error_max, abs(residuo));
            end
        end

        if error_max < tol
            break;
        end
    end
end

Ejemplos de aplicación

Caso 1: Ec. de Poisson con solución conocida

% Configuración del problema
xmin = 0; xmax = 1;
ymin = 0; ymax = 1;
nx = 51; ny = 51;

f = @(x,y) 2*pi^2 * sin(pi*x) .* sin(pi*y);
g = @(x,y) zeros(size(x));
sol_exacta = @(x,y) sin(pi*x) .* sin(pi*y);

[U_num, X, Y] = resolverEDP2D(f, g, xmin, xmax, ymin, ymax, nx, ny);
U_ref = sol_exacta(X, Y);

error_abs = abs(U_num - U_ref);
error_max = max(error_abs(:));

% Visualización
figure;
subplot(1,2,1);
surf(X, Y, U_num);
title('Solución numérica');

subplot(1,2,2);
surf(X, Y, error_abs);
title('Distribución del error');

Caso 2: Condiciones de frontera no homogéneas

xmin = -1; xmax = 1;
ymin = -1; ymax = 1;
nx = 41; ny = 41;

f = @(x,y) zeros(size(x));
g = @(x,y) x.^2 - y.^2;

U_sor = resolverSOR(f, g, xmin, xmax, ymin, ymax, nx, ny, 1.7, 500, 1e-6);

Análisis de convergencia

Se verifica el orden de covnergencia mediante refinamiento de malla:

tamaños = [11, 21, 41, 81];
errores = [];
pasos = [];

for n = tamaños
    h = 1/(n-1);
    [U, ~, ~] = resolverEDP2D(f, g, 0, 1, 0, 1, n, n);
    U_ex = sol_exacta(linspace(0,1,n), linspace(0,1,n));
    errores = [errores; max(abs(U(:) - U_ex(:)))];
    pasos = [pasos; h];
end

% Estimación del orden
for k = 1:length(tamaños)-1
    orden = log(errores(k)/errores(k+1)) / log(pasos(k)/pasos(k+1));
    fprintf('Orden estimado entre h=%.3f y h=%.3f: %.2f\n', pasos(k), pasos(k+1), orden);
end

Consideraciones técnicas

• Precisión: El esquema de cinco puntos alcanza precisión O(h²)
• Estabilidad: Los métodos directos son estables para sistemas moderados; los iterativos requieren ajuste del parámetro de relajación
• Eficiencia: La utilización de matrices dispersas reduce significativamente el uso de memoria