Newton Ralphson Iteration Getting Stuck on the While Loop

%%

%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.

Answer 1

[Alagu Sankar Esakkiappan]

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.