function [approximation, error] = newton_interpolation(x_values,
f_values, x_star)
% Calcula el polinomio interpolante de Newton de diferencias
divididas
% y devuelve la aproximación en x_star y el error de aproximación.
% Verifica si los tamaños de los vectores de valores de x y f son
iguales
if length(x_values) ~= length(f_values)
error('Los tamaños de los vectores de x y f deben ser
iguales.');
end
% Número de puntos de interpolación
n = length(x_values);
% Calcula las diferencias divididas
f_diff = f_values; % Inicializa con los valores de f
for j = 2:n
for i = n:-1:j
f_diff(i) = (f_diff(i) - f_diff(i-1)) / (x_values(i) -
x_values(i-j+1));
end
end
% Calcula la aproximación usando el polinomio interpolante de
Newton
approximation = f_diff(n);
for i = n-1:-1:1
approximation = approximation * (x_star - x_values(i)) +
f_diff(i);
end
% Calcula el error de aproximación
error = abs(f(x_star) - approximation);
end
Código en MATLAB
%Huaroc Loyola Luis
function [approximation, error] = simpsons_rule_approximation()
% Definir la función y los límites de integración
f = @(x) 1/sqrt(x^2 - 4);
a = 3;
b = 5;
% Número de subintervalos
n = 10;
% Calcular la longitud del subintervalo
delta_x = (b - a) / n;
% Calcular la aproximación usando la regla compuesta de Simpson
approximation = (delta_x / 3) * (f(a) + 4 * sum_odd_terms(f, a, n,
delta_x) + 2 * sum_even_terms(f, a, n, delta_x) + f(b));
% Calcular el error de aproximación
error = simpsons_rule_error(f, a, b, n);
% Mostrar resultados
disp('Aproximación:');
disp(approximation);
disp('Error de aproximación:');
disp(error);
end
function result = sum_odd_terms(f, a, n, delta_x)
result = 0;
for i = [Link]n-1
result = result + f(a + i * delta_x);
end
end
function result = sum_even_terms(f, a, n, delta_x)
result = 0;
for i = [Link]n-2
result = result + f(a + i * delta_x);
end
end
function error = simpsons_rule_error(f, a, b, n)
xi = linspace(a, b, n+1);
max_derivative = max(arrayfun(@(x) abs(fourth_derivative(f, x)),
xi));
error = ((b - a)^5) / (180 * n^4) * max_derivative;
end
function derivative = fourth_derivative(f, x)
h = 1e-5;
derivative = (f(x - 2*h) - 4*f(x - h) + 6*f(x) - 4*f(x + h) + f(x
+ 2*h)) / (h^4);
end
function [x] = gaussEliminacionPivoteoParcial(A, b)
[m, n] = size(A);
if m ~= n
error('La matriz A debe ser cuadrada.');
end
% Combinar matriz A y vector b
Ab = [A, b];
% Eliminación gaussiana con pivoteo parcial
for k = 1:n-1
% Encontrar el pivote máximo en la columna actual
[~, pivot_row] = max(abs(Ab(k:end, k)));
% Ajustar el índice de fila
pivot_row = pivot_row + k - 1;
% Intercambiar filas para tener el pivote en la posición
correcta
Ab([k, pivot_row], :) = Ab([pivot_row, k], :);
% Eliminación gaussiana
for i = k+1:n
factor = Ab(i, k) / Ab(k, k);
Ab(i, k:end) = Ab(i, k:end) - factor * Ab(k, k:end);
end
end
% Sustitución hacia atrás
x = zeros(n, 1);
for i = n:-1:1
x(i) = (Ab(i, end) - Ab(i, i+1:n) * x(i+1:n)) / Ab(i, i);
end
end
LUEGO EJECUTAMOS LA FUNCION PARA NUESTRO SISTEMA
function [t_values, y_values] = euler_modificado(tspan, y0, h,
tolerance)
% Método de Euler modificado para resolver ecuaciones
diferenciales ordinarias
% tspan: intervalo de tiempo [t_inicial, t_final]
% y0: condición inicial
% h: tamaño del paso
% tolerance: tolerancia para el error
% Inicialización de los vectores de tiempo y solución
t_values = tspan(1):h:tspan(2);
y_values = zeros(size(t_values));
y_values(1) = y0;
% Función diferencial
f = @(t, y) cos(2*t) + sin(3*t);
% Solución exacta para comparación
y_exact = @(t) (1/2) * sin(2*t) - (1/3) * cos(3*t) + (4/3);
% Iteración
for i = 1:length(t_values)-1
% Predicción usando Euler
y_pred = y_values(i) + h * f(t_values(i), y_values(i));
% Corrección usando el promedio de las derivadas
y_values(i+1) = y_values(i) + h * (f(t_values(i), y_values(i))
+ f(t_values(i+1), y_pred)) / 2;
% Calcular el error actual
error = abs(y_exact(t_values(i+1)) - y_values(i+1));
% Verificar la tolerancia
if error < tolerance
break;
end
end
end
% Parámetros del problema
tspan = [0, 1];
y0 = 1;
h = 0.1;
tolerance = 1e-5;
% Llamada al método de Euler modificado
[t_values, y_values] = euler_modificado(tspan, y0, h, tolerance);
% Muestra los resultados
plot(t_values, y_values, '-o', 'DisplayName', 'Aproximación');
hold on;
t_exact = linspace(tspan(1), tspan(2), 100);
y_exact = (1/2) * sin(2*t_exact) - (1/3) * cos(3*t_exact) + (4/3);
plot(t_exact, y_exact, '--', 'DisplayName', 'Solución exacta');
xlabel('t');
ylabel('y(t)');
title('Aproximación usando Euler modificado y solución exacta');
legend('Location', 'Best');
hold off;
