Why doesn't Newton's method calculate the correct temperature

function main
%UNTITLED Summary of this function goes here
%   Detailed explanation goes here
x0 = 80;
eps = 0.4;
c_s = 5.67*10^(-8);
s1 = 0.250;
s2 = 0.015;
lamda1 = 0.35;
lamda2 = 22.7;
Tw_1 = 1200;
T_l = 10;
Tw_1 = 1200;
T_l = 10;
lambda_L = 25.12;
betha = 3.543*10^(-3);
v = 144.0*10^(-7);
L = 7;
B = 15;              
A = B * L;           
g = 9.81;
Pr = 0.7095;
Ra = @(x)Pr*(g*L^3*betha*(x-T_l))/v^2;
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);
xsol = newton_raphson(func,x0)
alpha_k = alpha_k(xsol)
alpha_ges = alpha_k + (eps *c_s*(xsol.^4-T_l.^4)/(xsol-T_l))
Q_ges = (1/((s1/lamda1+s2/lamda2)+1/alpha_ges))*A*(Tw_1-T_l)
(1/(s1/lamda1+s2/lamda2))*(Tw_1-xsol) % heat conduction
alpha_k * (xsol-T_l) + (eps*c_s*xsol^4) % convection + thermal radiation

end
function xsol = newton_raphson(func,x0)
%UNTITLED Summary of this function goes here
%   Detailed explanation goes here
x(1) = x0;
maxiter = 200;
tol = 10^(-5);


for i = 1:maxiter
    
    diff_xi = (func(x(i)+1e-6)-func(x(i)))*1e6;
    if  diff_xi < tol
        fprintf('Pitfall hast occured a better initial guess\n');
        return;
        
    end
    
    x(i+1) = x(i) - func(x(i))/diff_xi;
    abs_error(i+1) = abs((x(i+1)-x(i))/x(i+1))*100;
    
    if abs(x(i+1) - x(i)) < tol && (1/(s1/lamda1+s2/lamda2))*(Tw_1-x(i+1)) == alpha_k(x(i+1)) * (x(i+1)-T_l) + (eps*c_s*x(i+1)^4) 
        fprintf('The Root has converged at x = %.10f\n', x(i+1));      
    else
        fprintf('Iteration no: %d,current guess x = %.10f, error = %.5f', i, x(i+1), abs_error(i+1));      
    end
    
end
xsol = x(end);
end

In the code above, the temperature, shown here as xsol, is to be calculated using the Newton method. Unfortunately, the heat transfer coefficient alpha_k also depends on the iterated xsol. I get a much too low temperature here. You can see that the code is not correct by the fact that heat conduction is not equal to heat convection+radiation.

(1/(s1/lamda1+s2/lamda2))*(Tw_1-xsol) % heat conduction
alpha_k * (xsol-T_l) + (eps*c_s*xsol^4) % convection + thermal radiation

Unfortunately, different values ​​come out here. I tried combining the if condition with an && to solve the problem. Unfortunately it does not work. I'm not very good at handling @functions at Matlab. I would be very grateful if someone could help me.

if abs(x(i+1) - x(i)) < tol && (1/(s1/lamda1+s2/lamda2))*(Tw_1-x(i+1)) == alpha_k(x(i+1)) * (x(i+1)-T_l) + (eps*c_s*x(i+1)^4) 

Answer 1

[Jon]

Is there a reason you are using your own solver here rather than the MATLAB fzero (for function of one variable) or fsolve for systems of non-linear equations?

Unless there is a specific reason not to use those, I would try using the built in solver. If you still have convergence issues, then the problem is most likely with how you have formulated the equations, and not the solver. You can then focus on what the problem is there.

I just skimmed your code briefly, but one aspect that might be good to look at more closely, is whether in fact your heat transfer coefficient is getting updated using the current value of your solution. You have to look carefully at how the anonymous functions are defined and the variables are scoped to make sure it isn't actually remaining constant, just using the initial guess, x0

[Aryo Aryanapour]

Unfortunately, I have to use a numerical method. I'm a student and it's mandatory. I have to be honest, I'm not studying computer science and I'm not very familiar with Matlab. I taught myself that little bit.

[Aryo Aryanapour]

I think that problem is due to dependency of alpha_k on xsoll. If I accept alpha_k as a constant value in the code, e.g. 4, then the correct solution comes out

[Aryo Aryanapour]

But I can't find it somehow in the code.

[Aryo Aryanapour]

It is a bit complicated. x is inserted into Ra and calculated, then alpha_k is calculated by Ra and at the end a new x value is determined using alpha_k. numerically using newoton. The new x-value must be reinserted into Ra and until Ra, alpha_k and x no longer change.

[Aryo Aryanapour]

Thank you anyway for your answer

[Jon]

You may have some other problems, but at the least I would suggest checking your function definitions for alpha_k and func. When I evaluate those using your code at your initial guess value of x0= 80 I get alpha_k(80) = 1.609e+20, and funct(80) = 1.4481e+22. Those are enoromous values, does that make physical sense?

[Aryo Aryanapour]

function main
%UNTITLED Summary of this function goes here
%   Detailed explanation goes here
x0 = 80;
eps = 0.4;
c_s = 5.67*10^(-8);
s1 = 0.250;
s2 = 0.015;
lamda1 = 0.35;
lamda2 = 22.7;
Tw_1 = 1200;
T_l = 10;
Tw_1 = 1200;
T_l = 10;
lambda_L = 25.12;
betha = 3.543*10^(-3);
v = 144.0*10^(-7);
L = 7;
B = 15;              
A = B * L;           
g = 9.81;
Pr = 0.7095;
Ra = @(x)Pr*(g*L^3*betha*(x-T_l))/v^2;
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);
xsol = newton_raphson(func,x0)
alpha_k = alpha_k(xsol)
alpha_ges = alpha_k + (eps *c_s*(xsol.^4-T_l.^4)/(xsol-T_l))
Q_ges = (1/((s1/lamda1+s2/lamda2)+1/alpha_ges))*A*(Tw_1-T_l)
(1/(s1/lamda1+s2/lamda2))*(Tw_1-xsol)
alpha_k * (xsol-T_l) + (eps*c_s*xsol^4)

end
function xsol = newton_raphson(func,x0)
%UNTITLED Summary of this function goes here
%   Detailed explanation goes here
x(1) = x0;
maxiter = 200;
tol = 10^(-5);


for i = 1:maxiter
    
    diff_xi = (func(x(i)+1e-6)-func(x(i)))*1e6;
    if  diff_xi < tol
        fprintf('Die Ableitung von x ist Null, ein anderes Startwert wählen\n');
        return;
        
    end
    
    x(i+1) = x(i) - func(x(i))/diff_xi;
    abs_error(i+1) = abs((x(i+1)-x(i))/x(i+1))*100;
    
    if abs(x(i+1) - x(i)) < tol 
        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);
end

so if I delete the superfluous && from the if condition, then I get the following values ​​with the code:

• xsol = 10 degrees Celsius (too little)
• alpha_k = 83.2230
• alpha_ges = 83.2231
• Q_ges = 1.7188e+05
• 1.6645e+03 % heat conduction
• 2.3027e-04 % convection + thermal radiation

[Jon]

It is hard to follow what your code is trying to do. I'm sorry I don't have time now to go through it and try to understand each step. Maybe someone else can assist you further.

[Jon]

The Newton-Raphson solves for a root (values that make function zero) of the function that you specify.

From your comments, it seems that you would like to find the temperature, x, that satisfies the energy balance:

conduction(x) = convection(x) + radiation(x)

or moving everything to the right hand side, you get the function you would like to find the root of

f(x) = -conduction(x) + convection(x) + radiation(x)

using the expressions for conduction, convection and radiation you have based upon the comments in your code you can define the function:

func = @(x) -(1/(s1/lamda1+s2/lamda2))*(Tw_1-x) + alpha_k(x) * (x-T_l) + (eps*c_s*x^4)

If I use this function, in place of yours, in your above code, I get:

Das Verfahren hat konvergiert bei x = 10.0000370032

Das Verfahren hat konvergiert bei x = 10.0000370032

xsol =

10.0000

alpha_k =

4.4982e+07

alpha_ges =

4.4982e+07

Q_ges =

1.7477e+05

ans =

1.6645e+03

ans =

1.6645e+03

So it seems like it converges ok and your checks are ok.

If you have reason to think the answer is not correct, you say you think it is too low, then you should probably recheck all of your expressions to make sure they are correct. Also, I would definitely be surprised if a celsius temperature should be used directly in an expression for the radiation. The radiation is typically in terms of Kelvin temperature to the 4th power, not Celsius as you have it.

Regarding your earlier termination criteria (original post)

if abs(x(i+1) - x(i)) < tol && (1/(s1/lamda1+s2/lamda2))*(Tw_1-x(i+1)) == alpha_k(x(i+1)) * (x(i+1)-T_l) + (eps*c_s*x(i+1)^4) 

This condition will probably never be met as it requires that the energy balance be achieved exactly (you use ==) Due to round off and limited precision representation of floating point numbers you should never look for exact equality, instead you should test that the two values are within some tolerance of each other.

A few other points:

You use create a variable named eps in your program, this should be avoided as eps is already a built in function which gives the machine precision.

You don't need to create a function main that takes no arguments and does not return anything. You can just delete off the first line of code

function main

and the last end statement, and run it as a MATLAB script. Just type the name of the file on the command line to run the script.