Implementación de Métodos Numéricos en MATLAB: Bisección, Newton, RKF45 y Tiro

Algoritmos de Métodos Numéricos en MATLAB

A continuación se detallan los códigos para la resolución de ecuaciones no lineales, integración adaptativa, ecuaciones diferenciales ordinarias mediante Runge-Kutta-Fehlberg y resolución de problemas de contorno por el método del tiro.

1. Método Híbrido: Newton-Raphson con Bisección

Este algoritmo combina la robustez del método de bisección con la rápida convergencia cuadrática del método de Newton-Raphson.

% Método Newton-Bisección
em = eps / 2;
a = 0;
b = 1;
fa = f(a);
fb = f(b);
x = (a + b) / 2;
fx = f(x);
[a, x, b], [fa, fx, fb]

a = x;
fa = fx;
dfx = df(x);
if abs(dfx) > em * abs(fx)
    x = x - fx / dfx;
    if a <= x & b >= x
        fx = f(x);
        [a, x, b], [fa, fx, fb]
    end
end

b = x;
fb = fx;
dfx = df(x);
if abs(dfx) > em * abs(fx)
    x = x - fx / dfx;
    if a <= x & b >= x
        fx = f(x);
        [a, x, b], [fa, fx, fb]
    end
end

2. Método Híbrido: Secante con Bisección

Aproximación de raíces combinando la estabilidad del intervalo de la bisección con la aproximación por pendientes del método de la secante.

% Método Secante-Bisección
em = eps / 2;
a = 1;
b = 2;
fa = f(a);
fb = f(b);
x0 = a;
x1 = b;
fx0 = fa;
fx1 = fb;
m = (fx1 - fx0) / (x1 - x0);
if abs(m) > em * abs(fx1)
    x = x1 - fx1 / m;
    if a <= x & b >= x
        fx = f(x);
        [a, x, b], [fa, fx, fb]
    end
end

b = x;
fb = fx;
x0 = x1;
x1 = x;
fx0 = fx1;
fx1 = fx;
m = (fx1 - fx0) / (x1 - x0);
if abs(m) > em * abs(fx1)
    x = x1 - fx1 / m;
    if a <= x & b >= x
        fx = f(x);
        [a, x, b], [fa, fx, fb]
    end
end

a = x;
fa = fx;
x0 = x1;
x1 = x;
fx0 = fx1;
fx1 = fx;
m = (fx1 - fx0) / (x1 - x0);
if abs(m) > em * abs(fx1)
    x = x1 - fx1 / m;
    if a <= x & b >= x
        fx = f(x);
        [a, x, b], [fa, fx, fb]
    end
end

3. Integración Numérica Adaptativa

Estrategia de cuadratura adaptativa que subdivide intervalos en función de los errores locales estimados.

% Código de Integración Numérica Adaptativa
x = [a, (a + b) / 2, b];
[i1, e1] = iehardy(x(1), x(2), f);
[i2, e2] = iehardy(x(2), x(3), f);
I = [i1, i2];
E = [e1, e2];
[s, k] = max(E);
itermax = 1000;

for iter = 1:itermax
    k1 = k + 1;
    xn = (x(k) + x(k1)) / 2;
    x = [x(1:k), xn, x(k1:end)];
    [i1, e1] = iehardy(x(k), xn, f);
    [i2, e2] = iehardy(xn, x(k1 + 1), f);
    I = [I(1:k - 1), i1, i2, I(k1:end)];
    E = [E(1:k - 1), e1, e2, E(k1:end)];
    error = sum(E);
    itg = sum(I);
    if error < ea
        [itg, error]
        return;
    else
        [s, k] = max(E);
    end
end

4. Método de Runge-Kutta-Fehlberg (RKF45)

Resolución de sistemas de ecuaciones diferenciales con control de paso adaptativo y estimación de error local.

% Método Runge-Kutta-Fehlberg (RKF45)
f = 'odesept2013';
t0 = 0;
y0 = [0; 1; 0];
tf = 30;
ea = 10^-12;
er = 10^-10;
[a, b, c, est] = rkf45a;
p1 = 1 / 5;
hmin = (tf - t0) / 10^6;
hn = max(hmin, (ea + er * norm(y0))^p1);
tn = t0;
yn = y0;
t = t0;
y = y0;

itmax = 1000;
for iter = 1:itmax
    [yn1, en] = rkf(f, tn, hn, yn, a, b, c, est);
    TOL = ea + er * norm(yn1);
    if en > TOL
        tn = tn + hn;
        if tn >= tf
            'Fin'
            return;
        end
        yn = yn1;
        t = [t, tn];
        y = [y, yn];
        rn = 0.9 * (TOL / en)^p1;
        if rn > 1.1 | rn < 1
            hn = min(5, rn) * hn;
        end
        if tn + hn > tf
            hn = tf - tn;
        end
    else
        [en, TOL]
        return;
    end
end

rn = 0.9 * (TOL / en)^p1;
hn = min(rn, 0.5) * hn;
[yn1, en] = rkf(f, tn, hn, yn, a, b, c, est);
TOL = ea + er * norm(yn1);
[en, TOL]

tn = tn + hn;
yn = yn1;
t = [t, tn];
y = [y, yn];
rn = 0.9 * (TOL / en)^p1;
hn = min(1, rn) * hn;
[yn1, en] = rkf(f, tn, hn, yn, a, b, c, est);
TOL = ea + er * norm(yn1);
[en, TOL]

5. Método del Tiro para Problemas de Contorno

Conversión de un problema de valores en la frontera en un problema de condiciones iniciales resuelto por iteraciones de la secante.

% Método de Tiro para Problemas de Contorno
em = eps / 2;
a = 1;
b = 2;
fa = tirof521(a);
fb = tirof521(b);
x0 = a;
x1 = b;
fx0 = fa;
fx1 = fb;
m = (fx1 - fx0) / (x1 - x0);
if abs(m) > em * abs(fx1)
    x = x1 - fx1 / m;
    [fx, tx, yx] = tirof521(x);
    [a, x, b], [fa, fx, fb]
end

x0 = x1;
x1 = x;
fx0 = fx1;
fx1 = fx;
a = x;
fa = fx;
m = (fx1 - fx0) / (x1 - x0);
if abs(m) > em * abs(fx1)
    x = x1 - fx1 / m;
    [a, x, b], [fx, tx, yx] = tirof521(x);
    [fa, fx, fb]
end

x0 = x1;
x1 = x;
fx0 = fx1;
fx1 = fx;
b = x;
fb = fx;
m = (fx1 - fx0) / (x1 - x0);
if abs(m) > em * abs(fx1)
    x = x1 - fx1 / m;
    [a, x, b], [fx, tx, yx] = tirof521(x);
    [fa, fx, fb]
end

x1 = x;
fx1 = fx;
b = x;
fb = fx;
m = (fx1 - fx0) / (x1 - x0);
if abs(m) > em * abs(fx1)
    x = x1 - fx1 / m;
    [a, x, b], [fx, tx, yx] = tirof521(x);
    [fa, fx, fb]
end