Implementaciones En Matlab

| |
IMPLEMENTACIONES EN MATLAB
* Cardano
function z = cardano (a)
if length(a)~= 4
error('Faltan (o sobran) coeficientes de la Ecuacion cubica')
z = NaN;
return
end
i = sqrt(-1);

a3 = a(1); a2 = a(2); a1 = a(3); a0 = a(4);

if a3 ~= 1
a0 = a0/a3; a1 = a1/a3; a2= a2/a3; a3 = 1;
end

t = a1 / 3 - a2 ^ 2 / 9; u = (a1 * a2 - 3 * a0) / 6 - a2 ^ 3 / 27;
d = t^3+ u^2; if abs(d) < 1e-16, d = 0; end

if (u+sqrt(d) < 0) & (d >= 0)
s1 = -abs((u+sqrt(d)))^(1/3);
else
s1 = (u + sqrt(d))^(1/3);
end
if (u-sqrt(d) < 0) & (d >= 0)
s2 = -abs((u-sqrt(d)))^(1/3);
else
s2= (u - sqrt(d))^(1/3);
end

if d > 0
disp('Una raiz real y dos complejas conjugadas')
elseif d < 0
disp('Todas las racies sonreales y diferentes')
else
disp('raices reales y almenos dos son igaules')
end

z(1) = (s1 + s2) - a2 /3;
z(2) = -(s1 + s2) / 2 - a2 / 3 + i * sqrt(3) / 2 * (s1 - s2);
z(3) = -(s1 + s2) / 2 - a2 / 3 - i * sqrt(3) / 2 * (s1 - s2);

z = sort(z);

function van = cubica(v)
van = v.^3-0.49836509*v.^2+0.0547*v-0.00166835;

function red= cubica(v)red=v.^3-0.46786509*v.^2-0.008017165*v-4.9083e-5

function soa= cubica(v)
soa=v.^3-0.46786509*v.^2+0.052982833*v-0.001319671

function pen= cubica(v)
pen=v.^3-0.448882377*v.^2+0.046950177*v-0.001073514

Ejecución

* Vander Waals
>> z=cardano([1,-0.49836509,0.0547,-0.00166835])'
Todas las racies son reales y diferentes
z =
0.0551
0.0843
0.3589

* Redlich – Kwong
>>z=cardano([1,-0.46786509,-0.008017165,-4.9083e-5])'
Una raiz real y dos complejas conjugadas
z =
-0.0084 + 0.0056i
-0.0084 - 0.0056i
0.4846

* Soave - Redlich – Kwong
>> z=cardano([1,-0.46786509,0.052982833,-0.001319671])
Todas las racies son reales y diferentes
z =
0.0348
0.1218
0.3113

* Peng – Robinson

>> z=cardano([1,-0.448882377,0.046950177,-0.001073514])'
Todas las racies son reales y diferentes
z =
0.0319
0.1094
0.3076

* Falsa Posición

function root= regfal (f, x1, x2)

epsilon = 1e-5;
Nmax= 50;

y1 = feval(f, x1); y2 = feval(f, x2);x3 = x2 -((x2 - x1)/(y2 - y1))* y2; y3 = feval(f, x3);

head;
headUp;

iter = 1; convergio = (1==0);
while (iter <= Nmax) & (~convergio)fprintf('%13d%13.6g%13.6g%13.6g%13.6g%13.6g%13.6g\n', iter, x1, x2, y1, y2, x3, y3)

convergio = abs(y3) < epsilon;

if ~convergio
iter = iter +1;
if y3 * y1 > 0
x1 = x3; y1 = y3;
else
x2 = x3; y2 = y3;
end
x3 = x2 -((x2 - x1)/(y2 - y1))* y2; y3 = feval(f, x3);
end
end

headDown;

ifconvergio
root = x3;
fprintf('\n%sValor de la raiz = %10.6g en %g iteraciones\n', b(7), root, iter)
else
root=NaN
fprintf('El metodo no convergio despues de %g iteraciones\n', Nmax)
end

function ans = b(n)
ans= blanks(n);

function ans = r (n)
ans = []; for i = 1:n, ans = [ans, '-']; end

function head
fprintf('\n%sAplicacion del METODO DE REGULA-FALSI\n\n',b(10))

function headUp
disp([b(9), 'iter', b(9), 'x1', b(9), 'x2', b(11), 'y1', b(12), 'y2', b(10), 'x3', b(12), 'y3'])
headDown

function headDown
fprintf('%s%s %s %s %s %s %s %s\n', b(9), r(6), r(12), r(12), r(12), r(12), r(12), r(12))

Ejecución

* Vander Waals
>> root= regfal('van', 0.6,0.4)

Aplicacion del METODO DE REGULA-FALSI

iter x1x2 y1 y2 x3 y3
------ ------------ ------------ ------------ ------------ ------------ ------------
1 0.6 0.4 0.0677402 0.00447324 0.385859 0.00268744
2 0.385859 0.4...