Newton Ralphson Iteration Getting Stuck on the While Loop

1 view (last 30 days)
Zaria Booth
Zaria Booth on 13 Dec 2021
%%
%Part 2 HW 8
clc; close all; clear all;
%stating the variables
MC = 8;
MH = 18;
Tr = 650; %reactant temperature
Y_Ycc = [0.70:0.1:1.40];
Ycc = MC + (MH/4);
Ymin = Ycc - (MC/2);
Y = (Y_Ycc)*(Ycc);
M_fuel = 114;
M_air = 28.97;
Rbar = 8.31451; %(kJ/kmol*K)
%new givens
P1=100; %(kPa) Pa=101.3
T1=300; %(K) Ta=288.2
rv=9.0; %compression ratio =V1/V2=V4/V3
V1=0.0005; %(m^3) ir 0.5L 6 cylinders=3.0 liters
eng_speed = 3000; %(RPM)
%find A,B,C for reactants and the N vaules for sum of hI for reactants
%for n=1:3
n=1:8;
% %new given
% y_fuel = (1)/(1+(4.76*Y(n)));
% y_air = (4.76*Y(n))/(1+(4.76*Y(n)));
% ans = y_fuel + y_air;
y_fuel = zeros(size(n));
y_air = zeros(size(n));
for p=1:numel(n)
y_fuel(p) = (1)/(1+(4.76*Y(n(p))));
y_air(p) = (4.76*Y(n(p)))/(1+(4.76*Y(n(p))));
end
%solving for n for reactants
n_tot = (P1*V1)/(Rbar*T1);
n_fuel = y_fuel*n_tot;
n_air = y_air*n_tot;
M_fuel_fix = n_fuel*M_fuel;
M_air_fix = n_air*M_air;
f_r = M_fuel_fix/M_air_fix;
%for Cpo
a_air = 27.5;
a_fuel = 38.4;
b_air = 0.0057;
b_fuel = 0.429;
%trying to find Tp by using newton's iteration
T2 = 5999; %initial guess is the highest I can go within 400<Tp<6000
iter = 0; %counts to number of iterations that take place
n=1; %declaring variable
%create a loop for the iteration of Tp
while n<=8
phi_air_T1 = ((a_air-Rbar)*(log(T1))) + (b_air*T1);
phi_fuel_T1 = ((a_fuel-Rbar)*(log(T1))) + (b_fuel*T1);
phi_air_T2 = ((a_air-Rbar)*(log(T2))) + (b_air*T2);
phi_fuel_T2 = ((a_fuel-Rbar)*(log(T2))) + (b_fuel*T2);
%find the part 3 of the equation for the products
m=1:8;
K1(m) = (-1*((y_air(n)*phi_air_T1)+(y_fuel(n)*phi_fuel_T1))) + (Rbar*log(1/rv));
FT2 = (y_air(n)*(phi_air_T2)) + (y_fuel(n)*phi_fuel_T2) + K1(n);
dphi_air_T2 = ((a_air-Rbar)/(T2)) + b_air;
dphi_fuel_T2 = ((a_fuel-Rbar)/(T2)) + b_fuel;
dFt2 = (y_air(n)*dphi_air_T2) + (y_fuel(n)*dphi_fuel_T2); %derivative of FT2
%sum the parts
sum_FT2 = sum(FT2,2);
sum_dFt2 = sum(dFt2,2);
%start the eqaution
T2_star = T2 - (FT2/dFt2);
tolerance = abs(T2_star-T2);
T2 = T2_star;
iter = iter + 1;
if tolerance < 1
T2_real(n)=T2_star;
%Z(n) = real(T2_real);
n=n+1; %allows the loop to do to n=13
end
end
%end
%for i=1:3
n=1:8;
P2 = zeros(size(n));
W = zeros(size(n));
for p=1:numel(n)
P2(p) = T2_real(n(p))*(P1/T1)*rv ;
part1(p) = n_air(n(p))* (((a_air-Rbar)*(T2_real(n(p))-T1)) + (0.5*b_air*((T2_real(n(p))^2)-(T1^2))));
part2(p) = n_fuel(n(p))* (((a_fuel-Rbar)*(T2_real(n(p))-T1)) + (0.5*b_fuel*((T2_real(n(p))^2)-(T1^2))));
W(p) = -1*((part1(n(p)))+(part2(n(p))));
end
answers = [Y/Ycc;T2_real;P2;W];
fprintf(' Y/Ycc T2 P2 W\n')
fprintf('--------------------------------------------------------------\n')
fprintf(' %1.5f | %5.7f | %5.7f | %5.7f\n',answers)
%begining of finding T3 and P3 sing T2_real
%ststing some givens
Hrp_To_Fuel = -5089100; %(kJ/kg) (kJ/kmol)=Hrp_To * M_fuel CHANGE BACK FOR HW
Urp_To = 281400;
Urp_To_Fuel = (Hrp_To_Fuel)-((2479.054)*((MH/4)-1));
%making the matrices for the A,B,and C values for the lower and upper limits
A = [299180 309070; %I=1 ==> CO
56835 93048; %I=2 ==> CO2
88923 154670; %I=3 ==> H2O
43388 127010; %I=4 ==> O2
31317 44639]; %I=5 ==> N2
B = [37.85 39.29; %I=1 ==> CO
66.27 68.58; %I=2 ==> CO2
49.36 60.43; %I=3 ==> H2O
42.27 46.25; %I=4 ==> O2
37.46 39.32]; %I=5 ==> N2
C = [-4571.9 -6201.9; %I=1 ==> CO
-11634.0 -16979.0; %I=2 ==> CO2
-7940.8 -19212.0; %I=3 ==> H2O
-6635.4 -18798.0; %I=4 ==> O2
-4559.3 -6753.4]; %I=5 ==> N2
D = [-31.10 -42.77; %I=1 ==> CO
-200.0 -220.4; %I=2 ==> CO2
-117.0 -204.6; %I=3 ==> H2O
-55.15 -92.15; %I=4 ==> O2
-34.82 -50.24]; %I=5 ==> N2
for n=1:8
n1 = 2*(Ycc-Y(n)); %n values for whe Y > Ycc and Ymin<Y<Ycc
n2 = 2*(Y(n)-Ymin);
n3 = MH/2;
n4 = 0;
n5 = 3.76*(Y(n));
n6 = 0;
n7 = MC;
n8 = MH/2;
n9 = Y(n)-Ycc;
n10 = 3.76*Y(n);
if (400<T2_real(n)) && (T2_real(n)<1600) %picks the correct A,B,C for reactants
Ar=A(:,1);
Br=B(:,1);
Cr=C(:,1);
elseif (1600<T2_real(n)) && (T2_real(n)<6000)
Ar=A(:,2);
Br=B(:,2);
Cr=C(:,2);
end
if (Y(n)<=Ycc) && (Ymin<=Y(n)) % picks N
N = [n1 n2 n3 n4 n5];
elseif (Y(n)>=Ycc)
N = [n6 n7 n8 n9 n10];
end
%part 2 of equation
for m=1:5
hI_1(n,m) = N(m)*(Ar(m) + (Br(m)*T2_real(n)) + (Cr(m)*log(T2_real(n))));
part3(n) = Urp_To_Fuel + (N(1)*Urp_To); %part 1 of the equation
end
%sum the values in reactants
sum_hI_1 = sum(hI_1,2);
part3 = part3'; %allows the matrixes to add/subtract
end
T3 = 5999; %initial guess is the highest I can go within 400<Tp<6000
iter = 0; %counts to number of iterations that take place
n=1; %declaring variable
%create a loop for the iteration of Tp
while n<=8
%stating my n values
n1 = 2*(Ycc-Y(n));
n2 = 2*(Y(n)-Ymin);
n3 = MH/2;
n4 = 0;
n5 = 3.76*(Y(n));
n6 = 0;
n7 = MC;
n8 = MH/2;
n9 = Y(n)-Ycc;
n10 = 3.76*Y(n);
if (400<T3) && (T3<1600)
Ap=A(:,1);
Bp=B(:,1);
Cp=C(:,1);
elseif (1600<T3) && (T3<6000)
Ap=A(:,2);
Bp=B(:,2);
Cp=C(:,2);
end
if (Y(n)<=Ycc) && (Ymin<=Y(n)) %picks N
N = [n1 n2 n3 n4 n5];
elseif (Y(n)>=Ycc)
N = [n6 n7 n8 n9 n10];
end
%find the part 3 of the equation for the products
for m=1:5
FT3(m) = N(m)*(Ap(m) + (Bp(m)*T3) + (Cp(m)*log(T3)));
dFt3(m) = N(m)*(Bp(m) + (Cp(m)/T3)); %derivative of FTp
part4(m,n) = N(m)*(Rbar*(T3-T2_real(n)));
end
%sum the parts
sum_FT3 = sum(FT3,2);
sum_dFt3 = sum(dFt3,2);
part4_in = part4';
sum_part4 = sum(part4_in,2);
%start the eqaution
T3_star = T3 - ((part3(n) + ((sum_FT3-sum_hI_1(n))-sum_part4(n)))/(sum_dFt3));
tolerance = abs(T3_star-T3);
T3 = T3_star;
iter = iter + 1;
if tolerance < 1
T3_real(n)=T3_star;
n=n+1; %allows the loop to do to n=13
end
end
for n=1:8
n1 = 2*(Ycc-Y(n)); %n values for whe Y > Ycc and Ymin<Y<Ycc
n2 = 2*(Y(n)-Ymin);
n3 = MH/2;
n4 = 0;
n5 = 3.76*(Y(n));
n6 = 0;
n7 = MC;
n8 = MH/2;
n9 = Y(n)-Ycc;
n10 = 3.76*Y(n);
if (Y(n)<=Ycc) && (Ymin<=Y(n)) % picks N
N = [n1 n2 n3 n4 n5];
elseif (Y(n)>=Ycc)
N = [n6 n7 n8 n9 n10];
end
sum_N(n,8) = (N(1)+N(2)+N(3)+N(4)+N(5))/(1+(4.76*Y(n)));
sum_N_sum = sum(sum_N,2);
sum_N_in = sum_N_sum';
%P3 = (P2(n)).*(sum_N_in(n)).*((T3_real(n))./(T2_real(n)));
y_co = (N(1))/(sum(N(1:5)));
y_co2 = (N(2))/(sum(N(1:5)));
y_h20 = (N(3))/(sum(N(1:5)));
y_o2 = (N(4))/(sum(N(1:5)));
y_n2 = (N(5))/(sum(N(1:5)));
end
n=1:8;
P3 = zeros(size(n));
for p = 1:numel(n)
P3(p) = (P2(n(p)))*(sum_N_in(n(p)))*((T3_real(n(p)))/(T2_real(n(p))));
end
for n=1:8
n1 = 2*(Ycc-Y(n)); %n values for whe Y > Ycc and Ymin<Y<Ycc
n2 = 2*(Y(n)-Ymin);
n3 = MH/2;
n4 = 0;
n5 = 3.76*(Y(n));
n6 = 0;
n7 = MC;
n8 = MH/2;
n9 = Y(n)-Ycc;
n10 = 3.76*Y(n);
if (Y(n)<=Ycc) && (Ymin<=Y(n)) % picks N
N = [n1 n2 n3 n4 n5];
elseif (Y(n)>=Ycc)
N = [n6 n7 n8 n9 n10];
end
sum_N(n,8) = (N(1)+N(2)+N(3)+N(4)+N(5))/(1+(4.76*Y(n)));
sum_N_sum = sum(sum_N,2);
sum_N_in = sum_N_sum';
%P3 = (P2(n)).*(sum_N_in(n)).*((T3_real(n))./(T2_real(n)));
y_co = (N(1))/(sum(N(1:5)));
y_co2 = (N(2))/(sum(N(1:5)));
y_h20 = (N(3))/(sum(N(1:5)));
y_o2 = (N(4))/(sum(N(1:5)));
y_n2 = (N(5))/(sum(N(1:5)));
n_tot_new = n_tot*((sum(N(1:5)))/(1+4.76*(Y(n))));
n_co = y_co*n_tot_new;
n_co2 = y_co2*n_tot_new;
n_h20 = y_h20*n_tot_new;
n_o2 = y_o2*n_tot_new;
n_n2 = y_n2*n_tot_new;
end
%creating the table
answers = [Y/Ycc;T3_real;P3];
fprintf('\n\n\n\n Y/Ycc T3 P3 \n')
fprintf('---------------------------------------------------\n')
fprintf(' %1.5f | %5.7f | %5.7f \n',answers)
for n=1:8
n1 = 2*(Ycc-Y(n)); %n values for whe Y > Ycc and Ymin<Y<Ycc
n2 = 2*(Y(n)-Ymin);
n3 = MH/2;
n4 = 0;
n5 = 3.76*(Y(n));
n6 = 0;
n7 = MC;
n8 = MH/2;
n9 = Y(n)-Ycc;
n10 = 3.76*Y(n);
if (400<T3_real(n)) && (T3_real(n)<1600) %picks the correct A,B,C for reactants
Ar=A(:,1);
Br=B(:,1);
Cr=C(:,1);
Dr=D(:,1)
elseif (1600<T3_real(n)) && (T3_real(n)<6000)
Ar=A(:,2);
Br=B(:,2);
Cr=C(:,2);
Dr=D(:,2);
end
if (Y(n)<=Ycc) && (Ymin<=Y(n)) % picks N
N = [n1 n2 n3 n4 n5];
elseif (Y(n)>=Ycc)
N = [n6 n7 n8 n9 n10];
end
%part 2 of equation
for m=1:5
phiT3(n,m) = (Br(m)*log(T3_real(n))) - ((Cr(m))/(T3_real(n))) + Dr(m);
y(m) = (N(m))/(sum(N(1:5)));
y_phiT3(n,m) = y(m)*(phiT3(n,m)); %part 1 of the equation
end
%sum the values in reactants
sum_phiT3 = sum(y_phiT3,2);
end
T4 = 5999; %initial guess is the highest I can go within 400<Tp<6000
iter = 0; %counts to number of iterations that take place
n=1; %declaring variable
%create a loop for the iteration of Tp
while n<=8
%stating my n values
n1 = 2*(Ycc-Y(n));
n2 = 2*(Y(n)-Ymin);
n3 = MH/2;
n4 = 0;
n5 = 3.76*(Y(n));
n6 = 0;
n7 = MC;
n8 = MH/2;
n9 = Y(n)-Ycc;
n10 = 3.76*Y(n);
if (400<T3) && (T3<1600)
Ap=A(:,1);
Bp=B(:,1);
Cp=C(:,1);
Dp=D(:,1);
elseif (1600<T3) && (T3<6000)
Ap=A(:,2);
Bp=B(:,2);
Cp=C(:,2);
Dp=D(:,2);
end
if (Y(n)<=Ycc) && (Ymin<=Y(n)) %picks N
N = [n1 n2 n3 n4 n5];
elseif (Y(n)>=Ycc)
N = [n6 n7 n8 n9 n10];
end
%find the part 3 of the equation for the products
for m=1:5
phiT4(m) = (Bp(m)*log(T4) - ((Cp(m))/(T4)) + Dp(m));
y_phiT4(m) = y(m)*(phiT4(m));
dFT4(m) = ((Bp(m))/(T4)) + ((Cp(m))/((T4)^2)); %derivative of FTp
y_dFT4(m) = y(m)*(dFT4(m));
end
%sum the parts
sum_phiT4 = sum(y_phiT4,2);
sum_dFT4 = sum(y_dFT4,2);
%equations
finalFT4 = (Rbar*log(rv)) - (Rbar*(log(T4)-log(T3_real(n)))) + (sum_phiT4(n)-sum_phiT3(n));
finaldFT4 = (sum_dFT4(n)) - (Rbar/T4);
%start the eqaution
T4_star = T4 - ((finalFT4(n))/(finaldFT4(n)));
tolerance = abs(T4_star-T4);
T4 = T4_star;
iter = iter + 1;
if tolerance < 1
T4_real(n)=T4_star;
n=n+1; %allows the loop to do to n=13
end
end
The Newton Ralpson iteration works for T3 and T2, but for some reason it won't for T4 unless Y_Ycc =1.0. Can anyone help me with this? The equation used for T4 is F(T4) = Rbar*ln(rv) - Rbar*ln(T4/T3) + sum(y(m)*(phi(T4(m))-phi(T3(m))) = 0.

Answers (1)

Alagu Sankar Esakkiappan
Alagu Sankar Esakkiappan on 17 Dec 2021
Hi Zaria,
I see that you're trying to find solution of an equation using Newton Raphson method. Convergence for Newton Raphson method highly depends on the initial starting point and is less consistent for Complex Roots. Off the bat, the main reason why iteration is stuck in T4 is because tolerance never reaches below the specified value of 1. i.e, tolerance consistently increases with each iteration and is also a complex number for T4 and consecutively n is never incremented. You may try checking the values for both sum_phiT3 and T3_real, if they are converged to expected values before solving for T4 ( Currently, both of them have converged to Complex Numbers before T4 ). Possible ways to find out the exact problem are intoducing breakpoints and checking step by step if converged values for all variables meet expected values.

Tags

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by