Taller de Metodos numericos

Tarea 1. Inciso 3.

Acosta Dominichetis J1, Carmona Sierra J1, Martínez Castellón L1, Upegui Hernandez M 1,
Teheran Manjarrez S1, Tovar Garces B1

1
Estudiantes de Ingeniería química de quinto semestre, Facultad de Ingeniería, Universidad de
Cartagena
3. Escoja un método abierto y otro cerrado para que realice un programa en el lenguaje
aprendido por usted y lo corra con la solución de estos ejercicios. Enviar el programa y
pantallazos de la corrida de los ejercicios.

Programa (método de bisección) Resolución de ejercicios inciso A

A continuación, se muestra el siguiente código implementado para la resolución de los ejercicios
a través del método de bisección, Nota: las capturas del código siendo ejecutado podrán ser
encontradas al final de este.

clc;
clear all;
syms x;
fprintf('Este es el método de bisección usando la ecuación de Van der Walls\n')
a = input ('Ingresa el valor de la constante a de Van der Walls (Pa*m^6/mol^2): ');
b = input ('Ingresa el valor de la constante b de Van der Walls (m^3/mol): ');
T = input ('Ingrese el valor de la temperatura (K): ');
P = input ('Ingrese el valor de la presión (Pa): ');
R = 8.314; % J/mol*k ; kPa*m^3/(mol*K)
f(x) = (P+a/x^2)*(x-b)-R*T;
xa = input('Ingrese el valor del limite inferior de su intervalo: ');
xb = input('Ingrese el valor del limite superior de su intervalo: ');
cc = 0.001;
Nmax = input('Ingrese el número de iteraciones: ');

fprintf('Ahora ingrese los siguientes datos para el gráfico:\n');

li = xa;
ls = xb;
 c = input('¿Cuantos puntos quiere graficar?: ');
p = linspace(li, ls, c);
y = double(subs(f(x),x,p));
plot(p, y, 'b');
grid on;
k = input('ingresa titulo del la variable independiente\n', 's');
b = input('ingresa titulo del la variable dependiente\n', 's')
xlabel(k);
ylabel(b);
if (f(xa) ~= 0 && f(xb) ~= 0)
if (f(xa) * f(xb) < 0)
i = 1;
while i <= Nmax

xr = vpa((xa + xb) / 2, 10);

n_i(i,1) = i;
n_xa(i,1) = xa;
n_xb(i,1) = xb;
n_xr(i,1) = double(xr);

if abs(f(xr)) <= cc
disp([num2str(double(xr)),' es la raíz']);
break;
else

if (f(xa) * f(xr) < 0)

xb = xr;
else
xa = xr;
 end
end
i = i + 1;
end
if i > Nmax
disp('Se alcanzó el número máximo de iteraciones sin converger');
end
else
disp(['No existe la raíz en el intervalo [', num2str(xa), ',', num2str(xb), '], elige un nuevo
intervalo']);
end
else
if f(xa) == 0
disp([num2str(xa),' es la raíz']);
else
disp([num2str(xb),' es la raíz']);
end
end
T = table(n_i, n_xa, n_xb, n_xr);
disp(T);
 - Para el Helio con los datos de P= 200000 Pa, T= 500 K, a = 0.0340 (Pa*m^6 / mol^2 ), b =
0.0234 * 10^-6

Para el 𝑯𝟐 con los datos P =250000 Pa, T = 450 K, a=0.244 (Pa*m^6 / mol^2 ), b = 0.0266 * 10^-
6 (m^3/mol)

Para el 𝑶𝟐 a P = 300000 Pa, T = 400 K , a = 1.36 (Pa*m^6 / mol^2 ) , b = 0.0318 * 10^-6 (m^3/mol)
Programa (Método de Newton Raphson) inciso B
A continuación se muestra el código empleado para la resolución del inciso B con el método
Newton Raphson, nota: Al final del mismo se encuentra evidencia de la ejecución del mismo

clc;
clear all;
syms x;
% Ingresar parámetros de ecuación virial
fprintf('Este es el método de Newton Raphson empleando ecuación virial\n');
B = input('Ingresa el valor de la constante B de la ecuación virial (cm³/mol): ');
C = input('Ingresa el valor de la constante C de la ecuación virial (cm⁶/mol²): ');
T = input('Ingrese el valor de la temperatura (K): ');
P = input('Ingrese el valor de la presión (bar): ');
R = 83.15;
f = (R*T/P)*(1+(B/x)+(C/x^2)) - x;
 % Código del gráfico
fprintf('Grafico de la función\n');
li = input('Ingresa el limite inferior de tu gráfica\n¿Desde qué valor la quieres?: ');
ls = input('Ingresa el limite superior de tu gráfica\n¿Hasta qué valor llegará?: ');
c = input('¿Cuantos puntos quiere graficar?: ');
p = linspace(li, ls, c);
y = double(subs(f,x,p));
plot(p, y, 'b');
grid on;
k = input('ingresa titulo del la variable independiente\n', 's');
b = input('ingresa titulo del la variable dependiente\n', 's')
xlabel(k);
ylabel(b);
repetir = true;

while repetir
df = diff(f,x);
xo = input('Ingrese un valor inicial:\n');
n = input ('Ingrese el número de iteraciones:\n');
tol = input ('Ingrese la tolerancia:\n');
for k = 1:n
fx = double(subs(f, x, xo));
dfx = double(subs(df, x, xo));

if dfx == 0
fprintf('la función no es iterable\n');
break;
end
 x1 = xo - fx/dfx;
fprintf('Iteración #%d: %.5f\n', k, x1);
if abs(x1-xo) < tol
fprintf('x%d es una raíz de la función : %.5f\n', k, x1);
break;
end
xo = x1;
end
respuesta = input('¿Es correcta la iteración? (si o no)', 's');
if strcmp(respuesta, 'si')
repetir = false;
end
z = P*x1/(R*T);
fprintf('Valor del factor de compresibilidad z: %.2f', z);
end
Para Isopropanol a T = 473.15 K, P = 10 bar, B = -388 (cm^3/mol) , C = -
26000 (cm^6/mol^2)
Programa (Método de Newton Raphson) inciso A
clc;
clear all;
% Solicitar los parámetros al usuario
a = input('Ingrese el valor de a (Pa·m^6/mol^2): ');
b = input('Ingrese el valor de b (m^3/mol): ');
P = input('Ingrese el valor de la presión P (Pa): ');
T = input('Ingrese el valor de la temperatura T (K): ');
R = 8.314;

% Limites de la grafica
v_min = input('Introduce el valor mínimo del volumen específico (m^3/mol) para la gráfica: ');
v_max = input('Introduce el valor máximo del volumen específico (m^3/mol) para la gráfica: ');

% f(v)
f = @(v) (P + a/v^2)*(v - b) - R*T;

% f '(v)
df = @(v) -2*a*(v - b)/v^3 + (P + a/v^2);

% Valores iniciales y tolerancia

v0 = input('Ingrese el valor inicial: ');
tol = input('Ingrese el valor de la tolerancia: ');
max_iter = 100;

% Inicialización de variables para guardar los resultados

iter_data = [];

% Solucion utilizando Newton Raphson

for iter = 1:max_iter
f_v0 = f(v0);
df_v0 = df(v0);
v_new = v0 - f_v0/df_v0;
error = abs(v_new - v0);

iter_data = [iter_data; iter, v0, f_v0, df_v0, v_new, error];

 v0 = v_new;

if error < tol

break;
end
end

% Mostrar tabla
tabla_resultados = array2table(iter_data, ...
'VariableNames', {'Iteración', 'x0', 'f(x0)', 'f''(x0)', 'xn', 'Error'});
disp(tabla_resultados);

% Verificar número máximo de iteraciones

if iter == max_iter
fprintf('No se alcanzó la convergencia después de %d iteraciones.\n', max_iter);
else
fprintf('El volumen específico convergió a %.6e m^3/mol en %d iteraciones.\n', v_new, iter);
end

v_values = linspace(v_min, v_max, 1000); % Rango de valores de volumen para graficar

f_values = arrayfun(f, v_values); % Evaluar la función en cada punto

% Grafica
figure;
plot(v_values, f_values, 'b-', 'LineWidth', 2);
xlabel('Volumen específico (v) [m^3/mol]');
ylabel('f(v)');
title('Función de Van der Waals');
xlim([v_min v_max]); % Establecer los límites del eje x con los valores ingresados por el usuario
grid on;
-Para el Helio con los datos de P = 200000 Pa, T = 500K, a = 0.0340 Pa*m^6*mol^-2, b =
0.0234*10^-6 m^3*mol^-1.

-Para el Hidrogeno con los datos de P = 250000 Pa, T = 450K, a = 0.244 Pa*m^6*mol^-2, b =
0.0266*10^-6 m^3*mol^-1.
-Para el Oxigeno con los datos de P = 300000 Pa, T = 400K, a = 1.36 Pa*m^6*mol^-2, b =
0.0318*10^-6 m^3*mol^-1.
Metodo de biseccion punto 2

clc;
clear all;

% Solicitar los parámetros al usuario

P = input('Introduce el valor de la presión P (bar): ');
T_Celsius = input('Introduce la temperatura en grados Celsius: ');
T = T_Celsius + 273.15; % Convertir a Kelvin
R = 83.14; % Constante de los gases (bar*cm^3/(mol·K))
B = input('Introduce el valor de B (cm^3/mol): ');
C = input('Introduce el valor de C (cm^6/mol^2): ');

% Definir la función de Z = PV/RT = 1 + B/V + C/V^2, reorganizada para f(V)

f = @(V) (R*T/P)*(1+(B./V)+(C./V.^2)) - V;

% Solicitar los valores mínimos y máximos para el eje x (volumen específico)

v_min = input('Introduce el valor mínimo del volumen específico (cm^3/mol) para la gráfica: ');
v_max = input('Introduce el valor máximo del volumen específico (cm^3/mol) para la gráfica: ');

% Método de bisección
xa = v_min;
xb = v_max;
tol = input('Ingrese el valor de la tolerancia: ');
max_iter = 100;

% Inicialización de variables para guardar resultados

iter_data = [];
xr = 0;
error = inf;

% Verificar que el intervalo es válido

if f(xa) * f(xb) > 0
error('El intervalo no encierra una raíz. f(xa) * f(xb) debe ser menor que 0.');
end

for i = 1:max_iter
fa = f(xa); % Evaluar f(xa)
fb = f(xb); % Evaluar f(xb)

xr_old = xr; % Guardar xr anterior

xr = (xa + xb) / 2; % Calcular el nuevo xr (punto medio)
fr = f(xr); % Evaluar f(xr)

% Calcular el error
if i > 1
error = abs((xr - xr_old) / xr); % Error relativo
end

% Guardar los datos de la iteración

iter_data = [iter_data; i, xa, xb, fa, fb, fa * fb, xr, fr, fr * fa, error];

% Verificar si f(xr)*f(xa) es positivo o negativo

if fa * fr < 0
xb = xr; % f(xr) y f(xa) tienen signos opuestos
else
 xa = xr; % f(xr) y f(xa) tienen el mismo signo
end

% Verificar si el error es menor que la tolerancia

if error < tol
break;
end
end

% Mostrar los resultados en forma de tabla

tabla_resultados = array2table(iter_data, ...
'VariableNames', {'Iteración', 'xa', 'xb', 'f(xa)', 'f(xb)', 'f(xa)*f(xb)', 'xr', 'f(xr)', 'f(xr)*f(xa)', 'Error'});
disp(tabla_resultados);

% Verificar si se alcanzó el número máximo de iteraciones

if i == max_iter
fprintf('No se alcanzó la convergencia después de %d iteraciones.\n', max_iter);
else
fprintf('El volumen molar convergió a %.6e cm^3/mol en %d iteraciones.\n', xr, i);
end

% Calcular Z usando el volumen molar encontrado

Z = (P * xr) / (R * T);
fprintf('El factor de compresibilidad Z es: %.6f\n', Z);

% Código del gráfico

fprintf('Gráfico de la función\n');
c = input('¿Cuántos puntos quiere graficar?: ');
p = linspace(v_min, v_max, c);
y = f(p); % Evaluar la función en los puntos generados

% Graficar
figure;
plot(p, y, 'b');
grid on;
xlabel('V (cm^3/mol)');
ylabel('f(V)');
title('Gráfico de la función de Van der Waals');

Para Isopropanol a T = 200 °C, P = 10 bar, B = -388 cm^3*mol^-1, C = -26*10^-3

cm^6*mol^-2.