Newton iteration doesn't run smoothly

1 view (last 30 days)
宏宇
宏宇 on 7 Sep 2022
Commented: Matt J on 8 Sep 2022
The author hopes to obtain the correct temperature distribution through Newton iteration.
This program is composed of three equations, which looks more complicated, but it is not difficult to understand because of the large amount of calculation. I have been troubled by various bugs in this program for a long time, and I hope to get your help
function1
% Set T2 to x.
function F=funfqtotceshi7(x)
global R q_tot fq_tot y i A B C TM K E S m
T1=300;
% T2=x;
TH=300;
TC=20;
g=5.675E-8;
c=0.04;
q_rad=(T1^4-x^4)*g/(1/c+1/c-1);
C1P=(1.1666E-3)*(1E-3);
a=0.9;
q_gas=C1P*a*(T1-x);
Tm=(T1+x)/2;
k=0.017+(7E-6)*(800-Tm)+0.0228*log(Tm);
C2f=0.008*0.02;
D=1/3000;
q_sol=(C2f*k/D)*(T1-x);
q_tot=q_rad+q_gas+q_sol;
% the first round of 'q_tot'
%-----------------------------
% Solve the first round of temperature distribution xn, assuming that there are only 8 layers.
syms xn1 xn2 xn3 xn4 xn5 xn6
eq1=solve(((1.1666E-3)*(1E-3)*0.9+0.008*0.02*(0.017+(7E-6)*(800-(x+xn1)/2)+0.0228*log((x+xn1)/2))/(1/3000))*(x-xn1)+(x^4-xn1^4)*5.675E-8/(1/0.04+1/0.04-1)-q_tot==0,xn1);
eq2=solve(((1.1666E-3)*(1E-3)*0.9+0.008*0.02*(0.017+(7E-6)*(800-(xn1+xn2)/2)+0.0228*log((xn1+xn2)/2))/(1/3000))*(xn1-xn2)+(xn1^4-xn2^4)*5.675E-8/(1/0.04+1/0.04-1)-q_tot==0,xn2);
eq3=solve(((1.1666E-3)*(1E-3)*0.9+0.008*0.02*(0.017+(7E-6)*(800-(xn2+xn3)/2)+0.0228*log((xn2+xn3)/2))/(1/3000))*(xn2-xn3)+(xn2^4-xn3^4)*5.675E-8/(1/0.04+1/0.04-1)-q_tot==0,xn3);
eq4=solve(((1.1666E-3)*(1E-3)*0.9+0.008*0.02*(0.017+(7E-6)*(800-(xn3+xn4)/2)+0.0228*log((xn3+xn4)/2))/(1/3000))*(xn3-xn4)+(xn3^4-xn4^4)*5.675E-8/(1/0.04+1/0.04-1)-q_tot==0,xn4);
eq5=solve(((1.1666E-3)*(1E-3)*0.9+0.008*0.02*(0.017+(7E-6)*(800-(xn4+xn5)/2)+0.0228*log((xn4+xn5)/2))/(1/3000))*(xn4-xn5)+(xn4^4-xn5^4)*5.675E-8/(1/0.04+1/0.04-1)-q_tot==0,xn5);
eq6=solve(((1.1666E-3)*(1E-3)*0.9+0.008*0.02*(0.017+(7E-6)*(800-(xn5+xn6)/2)+0.0228*log((xn5+xn6)/2))/(1/3000))*(xn5-xn6)+(xn5^4-xn6^4)*5.675E-8/(1/0.04+1/0.04-1)-q_tot==0,xn6);
A=[T1 x eq1 eq2 eq3 eq4 eq5 eq6];
%the first round of temperature
n=length(A);
B=sym(zeros(1,7));
C=sym(zeros(1,7));
TM=sym(zeros(1,7));
K=sym(zeros(1,7));
E=sym(zeros(1,7));
R=sym(zeros(1,7));
y=sym(zeros(1,8));
for i=sym(1:n-1)
B(i)=(A(i).^4-A(i+1).^4)*g/(1/c+1/c-1); %Radiant heat transfer
C(i)=C1P*a*(A(i)-A(i+1)); %Gas heat conduction
TM(i)=(A(i)+A(i+1))/2; %average temperature
K(i)=0.017+(7E-6)*(800-TM(i))+0.0228*log(TM(i)); %Thermal conductivity
E(i)=(C2f*K(i)/D)*(A(i)-A(i+1)); %Solid heat conduction
R(i)=1/((B(i)+C(i)+E(i))/(A(i)-A(i+1))); %Thermal resistance distribution
S=sum(R); %total thermal resistance
end
for m=sym(1:7)
y(1)=TC;
y(m+1)=TC+((sum(R(8-m:7)))/S)*(TH-TC);
% the second round of temperature y
fq_tot=(y(8)^4-y(7)^4)*g/(1/c+1/c-1)+C1P*a*(y(8)-y(7))+(C2f*(0.017+(7E-6)*(800-(y(8)+y(7))/2)+0.0228*log((y(8)+y(7))/2))/D)*(y(8)-y(7));
% the second round of fq_tot
end
F=q_tot-fq_tot;
end
function2
function dF=dfunfqtotceshi4(x)
global R q_tot fq_tot y i A B C TM K E S m
syms x
% syms xn1 xn2 xn3 xn4 xn5 xn6
T1=300;
% assumption T2
% T2=x;
TH=300;
TC=20;
g=5.675E-8;
c=0.04;
q_rad=(T1^4-x^4)*g/(1/c+1/c-1);
C1P=(1.1666E-3)*(1E-3);
a=0.9;
q_gas=C1P*a*(T1-x);
Tm=(T1+x)/2;
k=0.017+(7E-6)*(800-Tm)+0.0228*log(Tm);
C2f=0.008*0.02;
D=1/3000;
q_sol=(C2f*k/D)*(T1-x);
q_tot=q_rad+q_gas+q_sol;
syms xn1 xn2 xn3 xn4 xn5 xn6
% the first round of 'q_tot'
%-----------------------------
% Solve the first round of temperature distribution xn, assuming that there are only 8 layers.
eq1=solve(((1.1666E-3)*(1E-3)*0.9+0.008*0.02*(0.017+(7E-6)*(800-(x+xn1)/2)+0.0228*log((x+xn1)/2))/(1/3000))*(x-xn1)+(x^4-xn1^4)*5.675E-8/(1/0.04+1/0.04-1)-q_tot==0,xn1);
eq2=solve(((1.1666E-3)*(1E-3)*0.9+0.008*0.02*(0.017+(7E-6)*(800-(xn1+xn2)/2)+0.0228*log((xn1+xn2)/2))/(1/3000))*(xn1-xn2)+(xn1^4-xn2^4)*5.675E-8/(1/0.04+1/0.04-1)-q_tot==0,xn2);
eq3=solve(((1.1666E-3)*(1E-3)*0.9+0.008*0.02*(0.017+(7E-6)*(800-(xn2+xn3)/2)+0.0228*log((xn2+xn3)/2))/(1/3000))*(xn2-xn3)+(xn2^4-xn3^4)*5.675E-8/(1/0.04+1/0.04-1)-q_tot==0,xn3);
eq4=solve(((1.1666E-3)*(1E-3)*0.9+0.008*0.02*(0.017+(7E-6)*(800-(xn3+xn4)/2)+0.0228*log((xn3+xn4)/2))/(1/3000))*(xn3-xn4)+(xn3^4-xn4^4)*5.675E-8/(1/0.04+1/0.04-1)-q_tot==0,xn4);
eq5=solve(((1.1666E-3)*(1E-3)*0.9+0.008*0.02*(0.017+(7E-6)*(800-(xn4+xn5)/2)+0.0228*log((xn4+xn5)/2))/(1/3000))*(xn4-xn5)+(xn4^4-xn5^4)*5.675E-8/(1/0.04+1/0.04-1)-q_tot==0,xn5);
eq6=solve(((1.1666E-3)*(1E-3)*0.9+0.008*0.02*(0.017+(7E-6)*(800-(xn5+xn6)/2)+0.0228*log((xn5+xn6)/2))/(1/3000))*(xn5-xn6)+(xn5^4-xn6^4)*5.675E-8/(1/0.04+1/0.04-1)-q_tot==0,xn6);
A=[T1 x eq1 eq2 eq3 eq4 eq5 eq6];
%the first round of temperature
n=length(A);
B=sym(zeros(1,7));
C=sym(zeros(1,7));
TM=sym(zeros(1,7));
K=sym(zeros(1,7));
E=sym(zeros(1,7));
R=sym(zeros(1,7));
y=sym(zeros(1,8));
for i=sym(1:n-1)
B(i)=(A(i).^4-A(i+1).^4)*g/(1/c+1/c-1); %Radiant heat transfer
C(i)=C1P*a*(A(i)-A(i+1)); %Gas heat conduction
TM(i)=(A(i)+A(i+1))/2; %average temperature
K(i)=0.017+(7E-6)*(800-TM(i))+0.0228*log(TM(i)); %Thermal conductivity
E(i)=(C2f*K(i)/D)*(A(i)-A(i+1)); %Solid heat conduction
R(i)=1/((B(i)+C(i)+E(i))/(A(i)-A(i+1))); %Thermal resistance distribution
S=sum(R); %total thermal resistance
end
for m=sym(1:7)
y(1)=TC;
y(m+1)=TC+((sum(R(8-m:7)))/S)*(TH-TC);
% the second round of temperature y
fq_tot=(y(8)^4-y(7)^4)*g/(1/c+1/c-1)+C1P*a*(y(8)-y(7))+(C2f*(0.017+(7E-6)*(800-(y(8)+y(7))/2)+0.0228*log((y(8)+y(7))/2))/D)*(y(8)-y(7));
% the second round of fq_tot
end
F=q_tot-fq_tot;
dF=diff(F,x); %对F求导
end
Newton iteration
clear;clc
format long;
syms x0
x0=298; % Define the iteration initial value.
eps = 0.001; % Accuracy requirements
for i = 1:10000
F = double(subs(funfqtotceshi7(x0),'x' ,x0));
dF = double(subs(dfunfqtotceshi4(x0),'x' ,x0));
x = x0 - F/dF;
if (abs(x-x0)) ~= eps
break;
end
x0 = x; % update iteration result
end
disp('Solve results：');
x
disp('number of iterations：');
i
The author finds that x is around 279 by the band number method. But the Newton iteration only iterates once, and the obtained x is 206.
Many thanks in adavance,
Hopefully there is a solution for this problem.

Sign in to answer this question.

Answers (1)

Matt J
Matt J on 7 Sep 2022
Possibly, you meant to have,
if (abs(x-x0)) <= eps
  2 Comments
宏宇
宏宇 on 8 Sep 2022
Edited: 宏宇 on 8 Sep 2022
Thanks Matt
Referring to the changes you proposed, I have made changes to the code, but new errors appear.
equation 1
% Set T2 to x.
function F=funfqtotceshi7(x)
global R q_tot fq_tot y i A B C TM K E S m
T1=300;
% T2=x;
TH=300;
TC=20;
g=5.675E-8;
c=0.04;
q_rad=(T1^4-x^4)*g/(1/c+1/c-1);
C1P=(1.1666E-3)*(1E-3);
a=0.9;
q_gas=C1P*a*(T1-x);
Tm=(T1+x)/2;
k=0.017+(7E-6)*(800-Tm)+0.0228*log(Tm);
C2f=0.008*0.02;
D=1/3000;
q_sol=(C2f*k/D)*(T1-x);
q_tot=q_rad+q_gas+q_sol;
% the first round of 'q_tot'
%-----------------------------
% Solve the first round of temperature distribution xn, assuming that there are only 8 layers.
syms xn1 xn2 xn3 xn4 xn5 xn6
eq1=solve(((1.1666E-3)*(1E-3)*0.9+0.008*0.02*(0.017+(7E-6).*(800-(x+xn1)/2)+0.0228*log((x+xn1)/2))/(1/3000)).*(x-xn1)+(x.^4-xn1.^4)*5.675E-8/(1/0.04+1/0.04-1)-q_tot==0,xn1);
eq2=solve(((1.1666E-3)*(1E-3)*0.9+0.008*0.02*(0.017+(7E-6).*(800-(eq1+xn2)/2)+0.0228*log((eq1+xn2)/2))/(1/3000)).*(eq1-xn2)+(eq1.^4-xn2.^4)*5.675E-8/(1/0.04+1/0.04-1)-q_tot==0,xn2);
eq3=solve(((1.1666E-3)*(1E-3)*0.9+0.008*0.02*(0.017+(7E-6).*(800-(eq2+xn3)/2)+0.0228*log((eq2+xn3)/2))/(1/3000)).*(eq2-xn3)+(eq2.^4-xn3.^4)*5.675E-8/(1/0.04+1/0.04-1)-q_tot==0,xn3);
eq4=solve(((1.1666E-3)*(1E-3)*0.9+0.008*0.02*(0.017+(7E-6).*(800-(eq3+xn4)/2)+0.0228*log((eq3+xn4)/2))/(1/3000)).*(eq3-xn4)+(eq3.^4-xn4.^4)*5.675E-8/(1/0.04+1/0.04-1)-q_tot==0,xn4);
eq5=solve(((1.1666E-3)*(1E-3)*0.9+0.008*0.02*(0.017+(7E-6).*(800-(eq4+xn5)/2)+0.0228*log((eq4+xn5)/2))/(1/3000)).*(eq4-xn5)+(eq4.^4-xn5.^4)*5.675E-8/(1/0.04+1/0.04-1)-q_tot==0,xn5);
eq6=solve(((1.1666E-3)*(1E-3)*0.9+0.008*0.02*(0.017+(7E-6).*(800-(eq5+xn6)/2)+0.0228*log((eq5+xn6)/2))/(1/3000)).*(eq5-xn6)+(eq5.^4-xn6.^4)*5.675E-8/(1/0.04+1/0.04-1)-q_tot==0,xn6);
A=[T1 x eq1 eq2 eq3 eq4 eq5 eq6];
%the first round of temperature
n=length(A);
B=sym(zeros(1,n-1));
C=sym(zeros(1,n-1));
TM=sym(zeros(1,n-1));
K=sym(zeros(1,n-1));
E=sym(zeros(1,n-1));
R=sym(zeros(1,n-1));
y=sym(zeros(1,n));
for i=sym(1:n-1)
B(i)=(A(i).^4-A(i+1).^4)*g/(1/c+1/c-1); %Radiant heat transfer
C(i)=C1P*a*(A(i)-A(i+1)); %Gas heat conduction
TM(i)=(A(i)+A(i+1))/2; %average temperature
K(i)=0.017+(7E-6)*(800-TM(i))+0.0228*log(TM(i)); %Thermal conductivity
E(i)=(C2f*K(i)/D)*(A(i)-A(i+1)); %Solid heat conduction
R(i)=1/((B(i)+C(i)+E(i))/(A(i)-A(i+1))); %Thermal resistance distribution
S=sum(R); %total thermal resistance
end
for m=sym(1:7)
y(1)=TC;
y(m+1)=TC+((sum(R(8-m:7)))/S)*(TH-TC);
% the second round of temperature y
fq_tot=(y(8)^4-y(7)^4)*g/(1/c+1/c-1)+C1P*a*(y(8)-y(7))+(C2f*(0.017+(7E-6)*(800-(y(8)+y(7))/2)+0.0228*log((y(8)+y(7))/2))/D)*(y(8)-y(7));
% the second round of fq_tot
end
F=q_tot-fq_tot;
end
equation 2
function dF=dfunfqtotceshi4(x)
global R q_tot fq_tot y i A B C TM K E S m
syms x
% syms xn1 xn2 xn3 xn4 xn5 xn6
T1=300;
% assumption T2
% T2=x;
TH=300;
TC=20;
g=5.675E-8;
c=0.04;
q_rad=(T1^4-x^4)*g/(1/c+1/c-1);
C1P=(1.1666E-3)*(1E-3);
a=0.9;
q_gas=C1P*a*(T1-x);
Tm=(T1+x)/2;
k=0.017+(7E-6)*(800-Tm)+0.0228*log(Tm);
C2f=0.008*0.02;
D=1/3000;
q_sol=(C2f*k/D)*(T1-x);
q_tot=q_rad+q_gas+q_sol;
syms xn1 xn2 xn3 xn4 xn5 xn6
% the first round of 'q_tot'
%-----------------------------
% Solve the first round of temperature distribution xn, assuming that there are only 8 layers.
eq1=solve(((1.1666E-3)*(1E-3)*0.9+0.008*0.02*(0.017+(7E-6).*(800-(x+xn1)/2)+0.0228*log((x+xn1)/2))/(1/3000)).*(x-xn1)+(x.^4-xn1.^4)*5.675E-8/(1/0.04+1/0.04-1)-q_tot==0,xn1);
eq2=solve(((1.1666E-3)*(1E-3)*0.9+0.008*0.02*(0.017+(7E-6).*(800-(eq1+xn2)/2)+0.0228*log((eq1+xn2)/2))/(1/3000)).*(eq1-xn2)+(eq1.^4-xn2.^4)*5.675E-8/(1/0.04+1/0.04-1)-q_tot==0,xn2);
eq3=solve(((1.1666E-3)*(1E-3)*0.9+0.008*0.02*(0.017+(7E-6).*(800-(eq2+xn3)/2)+0.0228*log((eq2+xn3)/2))/(1/3000)).*(eq2-xn3)+(eq2.^4-xn3.^4)*5.675E-8/(1/0.04+1/0.04-1)-q_tot==0,xn3);
eq4=solve(((1.1666E-3)*(1E-3)*0.9+0.008*0.02*(0.017+(7E-6).*(800-(eq3+xn4)/2)+0.0228*log((eq3+xn4)/2))/(1/3000)).*(eq3-xn4)+(eq3.^4-xn4.^4)*5.675E-8/(1/0.04+1/0.04-1)-q_tot==0,xn4);
eq5=solve(((1.1666E-3)*(1E-3)*0.9+0.008*0.02*(0.017+(7E-6).*(800-(eq4+xn5)/2)+0.0228*log((eq4+xn5)/2))/(1/3000)).*(eq4-xn5)+(eq4.^4-xn5.^4)*5.675E-8/(1/0.04+1/0.04-1)-q_tot==0,xn5);
eq6=solve(((1.1666E-3)*(1E-3)*0.9+0.008*0.02*(0.017+(7E-6).*(800-(eq5+xn6)/2)+0.0228*log((eq5+xn6)/2))/(1/3000)).*(eq5-xn6)+(eq5.^4-xn6.^4)*5.675E-8/(1/0.04+1/0.04-1)-q_tot==0,xn6);
A=[T1 x eq1 eq2 eq3 eq4 eq5 eq6];
%the first round of temperature
n=length(A);
B=sym(zeros(1,n-1));
C=sym(zeros(1,n-1));
TM=sym(zeros(1,n-1));
K=sym(zeros(1,n-1));
E=sym(zeros(1,n-1));
R=sym(zeros(1,n-1));
y=sym(zeros(1,n));
for i=sym(1:n-1)
B(i)=(A(i).^4-A(i+1).^4)*g/(1/c+1/c-1); %Radiant heat transfer
C(i)=C1P*a*(A(i)-A(i+1)); %Gas heat conduction
TM(i)=(A(i)+A(i+1))/2; %average temperature
K(i)=0.017+(7E-6)*(800-TM(i))+0.0228*log(TM(i)); %Thermal conductivity
E(i)=(C2f*K(i)/D)*(A(i)-A(i+1)); %Solid heat conduction
R(i)=1/((B(i)+C(i)+E(i))/(A(i)-A(i+1))); %Thermal resistance distribution
S=sum(R); %total thermal resistance
end
for m=sym(1:7)
y(1)=TC;
y(m+1)=TC+((sum(R(8-m:7)))/S)*(TH-TC);
% the second round of temperature y
fq_tot=(y(8)^4-y(7)^4)*g/(1/c+1/c-1)+C1P*a*(y(8)-y(7))+(C2f*(0.017+(7E-6)*(800-(y(8)+y(7))/2)+0.0228*log((y(8)+y(7))/2))/D)*(y(8)-y(7));
% the second round of fq_tot
end
F=q_tot-fq_tot;
dF=diff(F,x); %diff to F
end
clear;clc
format long;
syms x0
x0=298; % Define the iteration initial value.
eps = 0.001; % Accuracy requirements
for i = 1:10000
F = double(subs(funfqtotceshi7(x0),'x' ,x0));
dF = double(subs(dfunfqtotceshi4(x0),'x' ,x0));
x = x0 - F/dF;
if (abs(x-x0)) <= eps
break;
end
x0 = x; % update iteration result
end
disp('Solve results：');
x
disp('number of iterations：');
i
error use sym/solve>getEqns
Input argument contains an empty equation or variable.
error sym/solve ( 226 line)
[eqns,vars,options] = getEqns(varargin{:});
error dfunfqtotceshi4 ( 35 line)
eq2=solve(((1.1666E-3)*(1E-3)*0.9+0.008*0.02*(0.017+(7E-6).*(800-(eq1+xn2)/2)+0.0228*log((eq1+xn2)/2))/(1/3000)).*(eq1-xn2)+(eq1.^4-xn2.^4)*5.675E-8/(1/0.04+1/0.04-1)-q_tot==0,xn2);
error newton3 ( 56 line)
dF = double(subs(dfunfqtotceshi4(x0),'x' ,x0));
Do you know what is causing this?
Matt J
Matt J on 8 Sep 2022
Possibly,
eq2=solve(((1.1666E-3)*(1E-3)*0.9+0.008*0.02*(0.017+(7E-6).*(800-(eq1+xn2)/2)+0.0228*log((eq1+xn2)/2))/(1/3000)).*(eq1-xn2)+(eq1.^4-xn2.^4)*5.675E-8/(1/0.04+1/0.04-1)-q_tot==0,xn2);
didn't have a solution and solve() therefore returned an empty [] result. You should check with the debugger.

Sign in to comment.

Categories

Find more on Calculus in Help Center and File Exchange

Tags

Products


Release

R2022a

Community Treasure Hunt

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

Start Hunting!