can someone help me pls with this code?

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

Unfortunately, I cannot calculate alpha_k and Ra in relation to x in this code. x is iterated using Newton's method and the result is then recalculated in alpha_k and Ra and x. The whole thing takes place until the x-value remains constant.

The result is correct if x no longer changes and these two equations have the same value

(1/(s1/lamda1+s2/lamda2))*(Tw_1-xsol) % heat conduction 

alpha_k * (xsol-T_l) + (eps*c_s*xsol^4) % convection + thermal radiation

Answer 1

[Saarthak Gupta]

Hi Aryo,

It seems like you are facing an issue while implementing the Newton-Raphson method correctly.

I suspect that you are getting inaccurate results for “alpha_k” and “Ra” where the Newton-Raphson method is converging. This might be due to the existence of multiple solutions. The function in question has multiple roots, specifically, 10 and -10. If you initialize "x0" with a negative number, for instance, -100, the Newton-Raphson method will converge to the root -10 after a sufficient number of iterations.

Refer to the following code:

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);

% calculate root of func using fzero
x_true = fzero(func,x0,optimset("Display","iter"));

 
Search for an interval around 80 containing a sign change:
 Func-count    a          f(a)             b          f(b)        Procedure
    1              80   1.44808e+22            80   1.44808e+22   initial interval
    3         77.7373   1.32188e+22       82.2627   1.58201e+22   search
    5            76.8   1.27182e+22          83.2    1.6398e+22   search
    7         75.4745   1.20319e+22       84.5255    1.7239e+22   search
    9            73.6   1.11038e+22          86.4   1.84764e+22   search
   11          70.949   9.87411e+21        89.051   2.03248e+22   search
   13            67.2   8.29396e+21          92.8   2.31423e+22   search
   15         61.8981   6.35876e+21       98.1019   2.75523e+22   search
   17            54.4   4.16874e+21         105.6   3.46918e+22   search
   19         43.7961   2.01761e+21       116.204   4.67419e+22   search
   21            28.8   4.50298e+20         131.2   6.81072e+22   search
   23         7.59227   3.34882e+18       152.408    1.0815e+23   search
   24           -22.4   -4.2743e+20       152.408    1.0815e+23   search
 
Search for a zero in the interval [-22.4, 152.408]:
 Func-count    x          f(x)             Procedure
   24           -22.4   -4.2743e+20        initial
   25        -21.7118  -3.86742e+20        interpolation
   26        -15.1941   -1.0826e+20        interpolation
   27        -15.1941   -1.0826e+20        bisection
   28        -14.1832  -8.03317e+19        interpolation
   29        -11.3014  -1.93896e+19        interpolation
   30        -11.3014  -1.93896e+19        bisection
   31        -9.62082   4.79324e+18        interpolation
   32        -9.95392   6.02432e+17        interpolation
   33        -10.0004  -4.75135e+15        interpolation
   34             -10    2.1968e+13        interpolation
   35             -10   7.94657e+08        interpolation
   36             -10             0        interpolation
 
Zero found in the interval [-22.4, 152.408]

% initialize x0 to a negative value
x = -100;
x_old = 1000;
iter = 0;
while abs(x_old-x) > 1e-10 && iter <= 20 % x ~= 0  
    x_old = x;
    dfunc = (func(x+1e-6)-func(x))*1e6;    % derivative of func
    x = x - func(x)/dfunc;
    iter = iter + 1;  
    fprintf('Iteration %d: x=%.18f, err=%.18f\n', iter, x, x_true-x);  
    pause(1);  
end 

Iteration 1: x=-65.862068404735197191, err=55.862068404735197191
Iteration 2: x=-43.270790174297857789, err=33.270790174297857789
Iteration 3: x=-28.477982878349731521, err=18.477982878349731521
Iteration 4: x=-19.052577706105886080, err=9.052577706105886080
Iteration 5: x=-13.475534197378028267, err=3.475534197378028267
Iteration 6: x=-10.793998287207900333, err=0.793998287207900333
Iteration 7: x=-10.056333899512672758, err=0.056333899512672758
Iteration 8: x=-10.000314686023203947, err=0.000314686023203947
Iteration 9: x=-10.000000009870499085, err=0.000000009870499085
Iteration 10: x=-9.999999999999998224, err=-0.000000000000001776
Iteration 11: x=-10.000000000000000000, err=0.000000000000000000

ak = alpha_k(x_true)

ak = 1.3134e+19

alpha_ges = ak + (eps *c_s*(x_true.^4-T_l.^4)/(x_true-T_l))

alpha_ges = 1.3134e+19

Q_ges = (1/((s1/lamda1+s2/lamda2)+1/alpha_ges))*A*(Tw_1-T_l)

Q_ges = 1.7477e+05

(1/(s1/lamda1+s2/lamda2))*(Tw_1-x_true) % heat conduction

ans = 1.6924e+03

ak * (x_true-T_l) + (eps*c_s*x_true^4) % convection + thermal radiation

ans = -2.6269e+20

I have used the "fzero" function rather than the Newton-Raphson method you implemented. "fzero" incorporates a mix of bisection, secant, and inverse quadratic interpolation techniques. To verify, consider comparing your outcomes with those obtained using this method.

Refer to the following MATLAB documentation for further reference:

• fzero: https://www.mathworks.com/help/matlab/ref/fzero.html

Hope this helps!

Best regards,

Saarthak