How can I find the last alpha_k value?

function main % Hauptfunktion
x0 = 200;       % Startwert für die erste Iteration. Ermittelt anhand des Plots der Funktion
eps = 0.4;      % Emissionsverhältnis 
c_s = 5.67*10^(-8);  % Strahlungskonstante
c12 = eps * c_s;     %
alpha_k = 4;         % Wärmeübergangskoeffizient bei ruhender Luft
s1 = 0.250;          % Wanddicke der feruerfesten Stein
s2 = 0.015;          % Wandddicke der äußeren Stahlwand
lamda1 = 0.35;       % Wärmeleitfähigkeit der feuerfesten Stein
lamda2 = 22.7;       % Wärmeleitfähigkeit der Stahlwand
Tw_1 = 1200;         % Wandtemperatur der feuerfesten Stein
T_l = 10;            % temperatur der Umgebungsluft
Pr = 0.7095;
lambda_L = 25.12;
betha = 3.543*10^(-3);
v = 144.0*10^(-7);
L = 7;
B = 15;              
A = B * L;           
g = 9.81;
Ra = @(x)Pr*(g*L^3*betha*(x-T_l))/v^2;
q1=0.825;
q2=0.387/((1+((0.492/Pr)^9/16))^(8/27));

q3=lambda_L/L;

alpha_k=@(x) ((q1+q2*(Ra(x)^(1/6)))^2)*q3;
%alpha_k = @(x)((0.825+(0.387*Ra(x)^1/6)/(1+(0.492/Pr)^9/16)^8/27)^2)*lambda_L/L;

func = @(x) eps * c_s * x^4 + (alpha_k(x) + 1/(s1/lamda1+s2/lamda2)) * x - ((1/(s1/lamda1+s2/lamda2)) * Tw_1 + alpha_k(x) * T_l); % Die Funktion
difft = @(x) (1713652349303737*x^3)/18889465931478580854784 - (138885763149071439966638298011785481*((6966792533989255*((2726811173664208678832570368*x)/66848762121410525 - 5453622347328417357665140736/13369752424282105)^(1/6))/18014398509481984 + 33/40))/(736008079697311707955200*((2726811173664208678832570368*x)/66848762121410525 - 5453622347328417357665140736/13369752424282105)^(5/6)) + (628*((6966792533989255*((2726811173664208678832570368*x)/66848762121410525 - 5453622347328417357665140736/13369752424282105)^(1/6))/18014398509481984 + 33/40)^2)/175 + (138885763149071439966638298011785481*x*((6966792533989255*((2726811173664208678832570368*x)/66848762121410525 - 5453622347328417357665140736/13369752424282105)^(1/6))/18014398509481984 + 33/40))/(7360080796973117079552000*((2726811173664208678832570368*x)/66848762121410525 - 5453622347328417357665140736/13369752424282105)^(5/6)) + 31780/22721;
xsol = newton_raphson(func, difft, x0)  
wl = (1/(s1/lamda1+s2/lamda2))*(Tw_1-xsol)          % Wärmestrom für Wärmeleitung
wsk = alpha_k(xsol) * (xsol-T_l) + (eps*c_s*xsol^4) % Wärmestrom für Konvektion + Strahlung
end

function xsol = newton_raphson(func, difft, x0)
x(1) = x0;
maxiter = 500;       % maximale Anzahl an Iterationen
Abbr = 10^(-5);      % Abbruchkriterium
eps = 0.4;           % Emissionsverhältnis 
c_s = 5.67*10^(-8);  % Strahlungskonstante
c12 = eps * c_s;     %
alpha_k = 4;         % Wärmeübergangskoeffizient bei ruhender Luft
s1 = 0.250;          % Wanddicke der feruerfesten Stein
s2 = 0.015;          % Wandddicke der äußeren Stahlwand
lamda1 = 0.35;       % Wärmeleitfähigkeit der feuerfesten Stein
lamda2 = 22.7;       % Wärmeleitfähigkeit der Stahlwand
Tw_1 = 1200;         % Wandtemperatur der feuerfesten Stein
T_l = 10;            % temperatur der Umgebungsluft
H = 7;               % Höhe
B = 15;              % Breite
A = B * H;           % Fläche

for i = 1:maxiter
    
    if  difft(x(i)) < Abbr % Ünerprüfung, ob die Ableitung von x Null ergibt. Falsches Startwert gewählt. Dadurch ergibt sich eine horizontale Tangente, welche die X-Achse nicht schneidet.Es kann keine weitere Nährungswert für die Nullstelle ausgeliefert weden.
        fprintf('Die Ableitung von x ist Null, ein anderes Startwert wählen\n');
        return;
        
    end
    
    
    x(i+1) = x(i) - func(x(i))/difft(x(i)); % Berechung mithilfe von Newton Verfahren (Hier befindet sich die Formel der Newton Verfahren)
    abs_error(i+1) = abs((x(i+1)-x(i))/x(i+1))*100;
    
    if abs(x(i+1) - x(i)) < Abbr
        fprintf('Das Verfahren hat konvergiert bei x = %.10f\n', x(i+1));      
    else
        fprintf('Iteration Nr: %d, aktuelle Schätzwert x = %.10f, error = %.5f, ', i, x(i+1), abs_error(i+1));      
    end
    
end
xsol = x(end);  % Oberflächentemperatur der Außenwand
alpha_k; % therefore I get the value 4. But the alpha k value changes after the iteration and is calculated depending on x. I need this value
alpha_k = alpha_k(xsol); % I need this value
alpha_ges = alpha_k + eps *c_s*(xsol^4-(T_l+273.15)^4)/(xsol-(T_l+273.15)) % Wärmeübergangskoeffizient-Gesamt
Q_ges = (1/((s1/lamda1+s2/lamda2)+1/alpha_ges))*A*(Tw_1-T_l) % Wärmestrom-Gesamt

end
Array indices must be positive integers or logical values.

Error in Copy_of_ruhender_Luft_rev1>newton_raphson (line 75)
alpha_k = alpha_k(xsol);

Error in Copy_of_ruhender_Luft_rev1 (line 32)
xsol = newton_raphson(func, difft, x0)
 

1) Can someone please help me to see the last value for alpha k (xsol) ? I need the alpha value calculated depending on x.

2) I would also like to put this code in a live script. Unfortunately it does not work. I want to keep the code but change some parameters and then put both codes (original and parameter modified ) in one script. can someone please tell me how to do this? Thanks in advance

[Torsten]

The value for xsol is passed to function main after successful Newton iteration.

Why do you want to calculate alpha_k twice : once in newton_raphson and once in main ?

But if this is necessary for some reason, I'd pass the function handle alpha_k to newton_raphson together with func, difft and x0, and remove the line alpha_k = 4 there.

Concerning the live script I have no experience - somebody else might help you in this respect.

[Aryo Aryanapour]

Hello Torsten,

thanks for the answer. I just want to know what the value of alpha_k is at the end of the iteration of x. Because alpha_k is calculated depending on x. If alpha_k ends up being 4, that's wrong. Because I overwrite the apha_k with the iteration. An X value is assumed, then an alpha is calculated with it, then a new x is calculated using Newton's method and then a new alpha_k with it.

[Aryo Aryanapour]

or am I misunderstanding something ?

[Torsten]

The function func has integrated the expression for alpha_k you defined in main as

alpha_k= @(x) ...

The value 4 for alpha_k set in newton_raphson is never used.

If you comment out this line

alpha_k = 4

you'll see that you arrive at the same result for xsol as before.

[Aryo Aryanapour]

Thank you so much Torsten. That's correct. The alpha_K = 4 was not used anywhere. This confused me a bit and I think this main function bothers me to run this code as a live script.

Answer 1

[Jon]

I was just going through your code and made a corrected version. I see now that Torsten has also given you a comment that echos what I have put in the attached file.

Your main problem is that you are reusing the same name alpha_k for a variable, which you assign the value of 4, and also as a function handle. I didn't actually see any place that you used the value of 4 assigned to alpha_k as a variable so I just commented that out.

I also added alpha_k as an argument to your function newton_raphson so that it would know what the definition of this function was.

Finally as I had told you in your earlier post on this same topic https://www.mathworks.com/matlabcentral/answers/1630740-why-doesn-t-newton-s-method-calculate-the-correct-temperature?s_tid=srchtitle you don't need to have the function main surrounding the start of the script. So I removed that. That might have been related to why this wasn't running as a live script. Since you defined the function main, you have to have a script or command line that calls main. You can't run a function with the run command. Anyhow I just copied and pasted my modified script into the live editor and it now runs fine. Also attached.

[Aryo Aryanapour]

hey joe

thank you very much. The code works and I can also run it in the live script.

[Jon]

Glad to hear this worked for you. One last thing, while you are cleaning up your code, as I mentioned previously you should use a different name for your variable that you call eps. There is already a MATLAB function eps which gives the machine precision. Renaming to a local variable can cause problems in the future if you actually try to use the MATLAB function, also it is confusing to others who may think you are invoking the MATLAB eps function.

[Aryo Aryanapour]

thank you very much Jon

yes That's right, you mentioned it before. I had so many variables. At some point I ran out of ideas. :)

I will define eps differently.

[Torsten]

What about "epsilon" ? Sounds nice :-)

[Aryo Aryanapour]

HAHA yes why not. I have also defined alpha as alpha.

[Jon]

Just don't use omicron :)

Answer 2

[Benjamin Thompson]

See the article "Anonymous Functions" in the MATLAB help. You must pass all parameters to the anonymous function as inputs, so maybe

Ra = @(x,Pr,g,L,betha,T_1,v) Pr*(g*L^3*betha*(x-T_l))/v^2;

Unless you want to declare all the parameters as globals.

Also, do not overwrite your function handle for alpha_k. Try:

alpha_k_output = alpha_k(xsol);

alpha_k will also need all those parameters that Ra needs unless you want them to be global. Try a smaller example to debug how you are using alpha_k and Ra first before trying to use them in the newton_raphson local function.

[Aryo Aryanapour]

Mr Thompson,

Many thanks for your help.