can someone help me pls with this code?

1 次查看（过去 30 天）
Aryo Aryanapour
Aryo Aryanapour 2022-1-17
编辑： Saarthak Gupta 2023-12-24
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

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

回答（1 个）

Saarthak Gupta
Saarthak Gupta 2023-12-24
编辑：Saarthak Gupta 2023-12-24
Ran in:
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:
Hope this helps!
Best regards,
Saarthak

类别

在 Help Center 和 File Exchange 中查找有关 Calculus 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by