Code does not update the array for each iteration of for loops

2 次查看（过去 30 天）
Matthew Palermo
Matthew Palermo 2023-11-26
评论： Matthew Palermo 2023-11-26
My code is not filling updating the temperature functino array (100x100) for each iteration of the for loops. Each column instead is a repeat of the next one. I understand that for the first couple of columns the answer is NaN becausre the solution should be -infinity, but I am not sure why it is simply repeating the previous results. Any help is much appreciated. I haave posted both functions and the script for clarity. Thanks!
function BB = feigR13(num,R,Bi)
x=R-1;
Cc=Bi*x;
p=1.5;
%Preallocation
BB = zeros(1,num);
for j= 1:num
R11=j*pi/x+0.155*(x)*(j-2)/j^2/(x+2.32)^1.8- 4*(j-1)*(0.062*x/(x+2)^1.56-0.0011)/j^2;
R12=(j-1/2)*pi/(x)- ((0.305+0.01*x^0.2+0.001*x)/(1+0.5*x)-0.12*(0.2*x^1)/(x^2+1.4)+0.067*x/(x^2+49))/j^1.5;
lambi=R11*Cc^p/(j*pi+Cc^p)+j*pi*R12/(j*pi+Cc^p);
lambi=round((lambi*10^10)/(10^10),-20);
x0 = lambi;
dx = pi/x/10;
for pp=1:6
UL = (x0 + dx);
LL = (x0 - dx);
ML = (UL + LL)/2;
f0=(-ML*bessely(1,ML*R)+Bi*bessely(0,ML*R))*besselj(0,ML)-(-ML*besselj(1,ML*R)+Bi*besselj(0,ML*R))*bessely(0,ML);
f1=(-LL*bessely(1,LL*R)+Bi*bessely(0,LL*R))*besselj(0,LL)-(-LL*besselj(1,LL*R)+Bi*besselj(0,LL*R))*bessely(0,LL);
f2=(-UL*bessely(1,UL*R)+Bi*bessely(0,UL*R))*besselj(0,UL)-(-UL*besselj(1,UL*R)+Bi*besselj(0,UL*R))*bessely(0,UL);
fp = (f2 - f1)/2/dx;
fpp = (f1 + f2 - 2*f0)/dx^2;
h = -f0/fp;
eps=-fpp*h^2/2/(fp + h*fpp);
x0 = x0 + h + eps;
dx = h^2;
if dx<10^-10
break
end
end
%Assign directly to output variable
BB(j)=x0;
end
end
clear all
close all
clc
%defining the dimensioned characteristics of the superheater element
A = 10; %systmm accuracy
k = 30; % BTU/hr-ft-F
h = 1000; % BTU/hr-ft^2-F
R1 = 0.137/2; %ft
R2 = 0.167/2; %ft
alpha = 0.0002; %ft^2/s thermal diffusivity for low-carbon steel (AISI 1010)
%function inputs
R = R2/R1;
rv = [0.01:0.01:1];
tv = [0.01:0.01:1];
Bi = (h*R1)/k;
r = length(rv);
t = length(tv);
TR13B01T0 = zeros(r,t);
for ir = 1:r
for it = 1:t
TR13B01T0(ir,it) = fDR13B01T0(rv(ir),tv(it),R,Bi,A);
end
end
plot(TR13B01T0(:,20),tv)
hold on
plot(TR13B01T0(:,50),tv)
hold on
plot(TR13B01T0(:,70),tv)
hold on
function [Td,qd]=fDR13B01T0(rv,tv,R,Bi,A)
srv=length(rv);
stv=length(tv);
Temp=zeros(stv,srv); % Preallocating arrays for speed
flux=zeros(stv,srv); % Preallocating arrays for speed
% calculate number of eigenvalues required to obtain solution with accuracy A
mmax1=floor((R-1)/pi*sqrt(A*log(10)/min(tv)));
bet=feigR13(mmax1,R,Bi);% Call the function to get eigenvalues
for ir = 1:srv % begin space loop
r=rv(ir);
term1=(Bi*R*log(r))/(1+Bi*R*log(R)); % steady-state temperature solution
term2=-(Bi*R/r)/(1+Bi*R*log(R)); % steady-state heat flux
%solution
for it=1:stv % begin time loop
t=tv(it);
mmax=floor((R-1)/pi*sqrt(A*log(10)/t));
term3 = 0;
term4 = 0;
for ii=1:mmax % begin summation loop
bt=bet(ii);
V0=-bt*besselj(1,bt*R)+Bi*besselj(0,bt*R);
S0=-bt*bessely(1,bt*R)+Bi*bessely(0,bt*R);
Nr = (S0*besselj(0,bt*r)-V0*bessely(0,bt*r));
L1 = (S0*besselj(1,bt*R)-V0*bessely(1,bt*R));
Denominator = (Bi^2+bt^2)*besselj(0,bt)*besselj(0,bt)-V0^2;
Norm = (Denominator)/(besselj(0,bt)*besselj(0,bt));
term3 = term3 - (pi*pi/2)*R*Bi*(exp(-bt*bt*t))*bt*Nr*L1/Norm;
Nm = (S0*besselj(1,bt*r)-V0*bessely(1,bt*r));
term4 = term4 - (pi*pi/2)*Bi*R*exp(-bt*bt*t)*bt*bt*Nm*L1/Norm;
end
T(ir) = abs(term1); % Steady-state temperature solution
Q(ir) = abs(term2); % Steady-state heat flux solution
Temp(it,ir) = term1 + term3; % Total temperature solution
flux(it,ir) = term2 + term4; % Total heat flux solution
end
end
Td=abs(Temp);
qd=abs(flux);
end
  2 个评论
Walter Roberson
Walter Roberson 2023-11-26
Ran in:
Re-arranging the code into more usable form:
%defining the dimensioned characteristics of the superheater element
A = 10; %systmm accuracy
k = 30; % BTU/hr-ft-F
h = 1000; % BTU/hr-ft^2-F
R1 = 0.137/2; %ft
R2 = 0.167/2; %ft
alpha = 0.0002; %ft^2/s thermal diffusivity for low-carbon steel (AISI 1010)
%function inputs
R = R2/R1;
rv = [0.01:0.01:1];
tv = [0.01:0.01:1];
Bi = (h*R1)/k;
r = length(rv);
t = length(tv);
TR13B01T0 = zeros(r,t);
for ir = 1:r
for it = 1:t
TR13B01T0(ir,it) = fDR13B01T0(rv(ir),tv(it),R,Bi,A);
end
end
plot(TR13B01T0(:,20),tv)
hold on
plot(TR13B01T0(:,50),tv)
hold on
plot(TR13B01T0(:,70),tv)
hold on
function BB = feigR13(num,R,Bi)
x=R-1;
Cc=Bi*x;
p=1.5;
%Preallocation
BB = zeros(1,num);
for j= 1:num
R11=j*pi/x+0.155*(x)*(j-2)/j^2/(x+2.32)^1.8- 4*(j-1)*(0.062*x/(x+2)^1.56-0.0011)/j^2;
R12=(j-1/2)*pi/(x)- ((0.305+0.01*x^0.2+0.001*x)/(1+0.5*x)-0.12*(0.2*x^1)/(x^2+1.4)+0.067*x/(x^2+49))/j^1.5;
lambi=R11*Cc^p/(j*pi+Cc^p)+j*pi*R12/(j*pi+Cc^p);
lambi=round((lambi*10^10)/(10^10),-20);
x0 = lambi;
dx = pi/x/10;
for pp=1:6
UL = (x0 + dx);
LL = (x0 - dx);
ML = (UL + LL)/2;
f0=(-ML*bessely(1,ML*R)+Bi*bessely(0,ML*R))*besselj(0,ML)-(-ML*besselj(1,ML*R)+Bi*besselj(0,ML*R))*bessely(0,ML);
f1=(-LL*bessely(1,LL*R)+Bi*bessely(0,LL*R))*besselj(0,LL)-(-LL*besselj(1,LL*R)+Bi*besselj(0,LL*R))*bessely(0,LL);
f2=(-UL*bessely(1,UL*R)+Bi*bessely(0,UL*R))*besselj(0,UL)-(-UL*besselj(1,UL*R)+Bi*besselj(0,UL*R))*bessely(0,UL);
fp = (f2 - f1)/2/dx;
fpp = (f1 + f2 - 2*f0)/dx^2;
h = -f0/fp;
eps=-fpp*h^2/2/(fp + h*fpp);
x0 = x0 + h + eps;
dx = h^2;
if dx<10^-10
break
end
end
%Assign directly to output variable
BB(j)=x0;
end
end
function [Td,qd]=fDR13B01T0(rv,tv,R,Bi,A)
srv=length(rv);
stv=length(tv);
Temp=zeros(stv,srv); % Preallocating arrays for speed
flux=zeros(stv,srv); % Preallocating arrays for speed
% calculate number of eigenvalues required to obtain solution with accuracy A
mmax1=floor((R-1)/pi*sqrt(A*log(10)/min(tv)));
bet=feigR13(mmax1,R,Bi);% Call the function to get eigenvalues
for ir = 1:srv % begin space loop
r=rv(ir);
term1=(Bi*R*log(r))/(1+Bi*R*log(R)); % steady-state temperature solution
term2=-(Bi*R/r)/(1+Bi*R*log(R)); % steady-state heat flux
%solution
for it=1:stv % begin time loop
t=tv(it);
mmax=floor((R-1)/pi*sqrt(A*log(10)/t));
term3 = 0;
term4 = 0;
for ii=1:mmax % begin summation loop
bt=bet(ii);
V0=-bt*besselj(1,bt*R)+Bi*besselj(0,bt*R);
S0=-bt*bessely(1,bt*R)+Bi*bessely(0,bt*R);
Nr = (S0*besselj(0,bt*r)-V0*bessely(0,bt*r));
L1 = (S0*besselj(1,bt*R)-V0*bessely(1,bt*R));
Denominator = (Bi^2+bt^2)*besselj(0,bt)*besselj(0,bt)-V0^2;
Norm = (Denominator)/(besselj(0,bt)*besselj(0,bt));
term3 = term3 - (pi*pi/2)*R*Bi*(exp(-bt*bt*t))*bt*Nr*L1/Norm;
Nm = (S0*besselj(1,bt*r)-V0*bessely(1,bt*r));
term4 = term4 - (pi*pi/2)*Bi*R*exp(-bt*bt*t)*bt*bt*Nm*L1/Norm;
end
T(ir) = abs(term1); % Steady-state temperature solution
Q(ir) = abs(term2); % Steady-state heat flux solution
Temp(it,ir) = term1 + term3; % Total temperature solution
flux(it,ir) = term2 + term4; % Total heat flux solution
end
end
Td=abs(Temp);
qd=abs(flux);
end
Walter Roberson
Walter Roberson 2023-11-26
Ran in:
Re-arranging the plotting as you had reversed x and y coordinates (but it does not make any difference to the overall problem)
%defining the dimensioned characteristics of the superheater element
A = 10; %systmm accuracy
k = 30; % BTU/hr-ft-F
h = 1000; % BTU/hr-ft^2-F
R1 = 0.137/2; %ft
R2 = 0.167/2; %ft
alpha = 0.0002; %ft^2/s thermal diffusivity for low-carbon steel (AISI 1010)
%function inputs
R = R2/R1;
rv = [0.01:0.01:1];
tv = [0.01:0.01:1];
Bi = (h*R1)/k;
r = length(rv);
t = length(tv);
TR13B01T0 = zeros(r,t);
for ir = 1:r
for it = 1:t
TR13B01T0(ir,it) = fDR13B01T0(rv(ir),tv(it),R,Bi,A);
end
end
plot(tv, TR13B01T0(:,20), 'displayname', string(tv(20)))
hold on
plot(tv, TR13B01T0(:,50), 'displayname', string(tv(50)))
hold on
plot(tv, TR13B01T0(:,70), 'displayname', string(tv(70)))
hold on
legend show
max((TR13B01T0(:,20) - TR13B01T0(:,50)).')
ans = 0
function BB = feigR13(num,R,Bi)
x=R-1;
Cc=Bi*x;
p=1.5;
%Preallocation
BB = zeros(1,num);
for j= 1:num
R11=j*pi/x+0.155*(x)*(j-2)/j^2/(x+2.32)^1.8- 4*(j-1)*(0.062*x/(x+2)^1.56-0.0011)/j^2;
R12=(j-1/2)*pi/(x)- ((0.305+0.01*x^0.2+0.001*x)/(1+0.5*x)-0.12*(0.2*x^1)/(x^2+1.4)+0.067*x/(x^2+49))/j^1.5;
lambi=R11*Cc^p/(j*pi+Cc^p)+j*pi*R12/(j*pi+Cc^p);
lambi=round((lambi*10^10)/(10^10),-20);
x0 = lambi;
dx = pi/x/10;
for pp=1:6
UL = (x0 + dx);
LL = (x0 - dx);
ML = (UL + LL)/2;
f0=(-ML*bessely(1,ML*R)+Bi*bessely(0,ML*R))*besselj(0,ML)-(-ML*besselj(1,ML*R)+Bi*besselj(0,ML*R))*bessely(0,ML);
f1=(-LL*bessely(1,LL*R)+Bi*bessely(0,LL*R))*besselj(0,LL)-(-LL*besselj(1,LL*R)+Bi*besselj(0,LL*R))*bessely(0,LL);
f2=(-UL*bessely(1,UL*R)+Bi*bessely(0,UL*R))*besselj(0,UL)-(-UL*besselj(1,UL*R)+Bi*besselj(0,UL*R))*bessely(0,UL);
fp = (f2 - f1)/2/dx;
fpp = (f1 + f2 - 2*f0)/dx^2;
h = -f0/fp;
eps=-fpp*h^2/2/(fp + h*fpp);
x0 = x0 + h + eps;
dx = h^2;
if dx<10^-10
break
end
end
%Assign directly to output variable
BB(j)=x0;
end
end
function [Td,qd]=fDR13B01T0(rv,tv,R,Bi,A)
srv=length(rv);
stv=length(tv);
Temp=zeros(stv,srv); % Preallocating arrays for speed
flux=zeros(stv,srv); % Preallocating arrays for speed
% calculate number of eigenvalues required to obtain solution with accuracy A
mmax1=floor((R-1)/pi*sqrt(A*log(10)/min(tv)));
bet=feigR13(mmax1,R,Bi);% Call the function to get eigenvalues
for ir = 1:srv % begin space loop
r=rv(ir);
term1=(Bi*R*log(r))/(1+Bi*R*log(R)); % steady-state temperature solution
term2=-(Bi*R/r)/(1+Bi*R*log(R)); % steady-state heat flux
%solution
for it=1:stv % begin time loop
t=tv(it);
mmax=floor((R-1)/pi*sqrt(A*log(10)/t));
term3 = 0;
term4 = 0;
for ii=1:mmax % begin summation loop
bt=bet(ii);
V0=-bt*besselj(1,bt*R)+Bi*besselj(0,bt*R);
S0=-bt*bessely(1,bt*R)+Bi*bessely(0,bt*R);
Nr = (S0*besselj(0,bt*r)-V0*bessely(0,bt*r));
L1 = (S0*besselj(1,bt*R)-V0*bessely(1,bt*R));
Denominator = (Bi^2+bt^2)*besselj(0,bt)*besselj(0,bt)-V0^2;
Norm = (Denominator)/(besselj(0,bt)*besselj(0,bt));
term3 = term3 - (pi*pi/2)*R*Bi*(exp(-bt*bt*t))*bt*Nr*L1/Norm;
Nm = (S0*besselj(1,bt*r)-V0*bessely(1,bt*r));
term4 = term4 - (pi*pi/2)*Bi*R*exp(-bt*bt*t)*bt*bt*Nm*L1/Norm;
end
T(ir) = abs(term1); % Steady-state temperature solution
Q(ir) = abs(term2); % Steady-state heat flux solution
Temp(it,ir) = term1 + term3; % Total temperature solution
flux(it,ir) = term2 + term4; % Total heat flux solution
end
end
Td=abs(Temp);
qd=abs(flux);
end

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

Walter Roberson
Walter Roberson 2023-11-26
lambi=round((lambi*10^10)/(10^10),-20);
Notice the round() covers all of (lambi*10^10)/(10^10) as its first parameter. But algebraically the multiplication and division cancel out, which would lead to round(lambi,-20) which is a request to round lambi to the nearest 10^20 . The existing code only differs from that to the extent that if lambi is sufficiently large then lambi*10^10 might not be exactly representable, so potentially the (lambi*10^10)/(10^10) changes the round-off of lambi.. or that lambi*10^10 might overflow to infinity. An example of a value that changes in going through *10^10/10^10 but is small enough to be affected by the round(-20), is lambi = 6.1479804705177653e+69
At the very least, that section needs to be heavily commented to describe why it is important to do that multiplying and diviging by 10^10 and that rounding to the nearest 10^20 .
The actual values that go into calculating lambi are relative constants -- the same for every pairs of r and t values. So before the odd rounding, lambi will always be roughly 7 2/3. There are values near that, such as 7.2664788935685998 that come out slightly different with the *10^10/10^10 process... but all values in that range round to 0 to the nearest 10^20 .
So... your lambi is rounding to 0. And then when you step through the UL, LL, ML calculation, ML will come out as 0. You have bessely(0,ML) but that is -inf ... so you have calculations involving infinities, and those will typically generate NaN.
The reason why you do not always generate NaN for all of your iterations, is that current value of t makes a difference in calculationg mmax1, and mmax1 is used as the upper limit of your for loop. But with t values beyond roughly 0.11 then mmax1 comes out as 0, so your loop is 1:0 which causes it not to be done at all, leaving the arrays still at the zero() that you initialized them to.
Basically your code is calculating whether t is more than some number greater than roughly 0.11, and if so then return 0, and otherwise return NaN.
  1 个评论
Matthew Palermo
Matthew Palermo 2023-11-26
Gotcha, I was overlooking that. I got this code from a scientific paper, but I have had issues using the authors' codes before.
Thanks!

请先登录，再进行评论。

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Particle & Nuclear Physics 的更多信息

标签

产品


版本

R2023b

Community Treasure Hunt

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

Start Hunting!

Translated by