Newton iteration doesn't run smoothly

Question

宏宇 2022년 9월 7일

0
링크

이 질문에 대한 바로 가기 링크

https://kr.mathworks.com/matlabcentral/answers/1798615-newton-iteration-doesn-t-run-smoothly

댓글: Matt J 2022년 9월 8일

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.

댓글 수: 0
이전 댓글 -2개 표시이전 댓글 -2개 숨기기

댓글을 달려면 로그인하십시오.

이 질문에 답변하려면 로그인하십시오.

Answer 1

Matt J 2022년 9월 7일

0
링크

이 답변에 대한 바로 가기 링크

https://kr.mathworks.com/matlabcentral/answers/1798615-newton-iteration-doesn-t-run-smoothly#answer_1046705

MATLAB Online에서 열기

Possibly, you meant to have,

if (abs(x-x0)) <= eps

댓글 수: 2
없음 표시없음 숨기기

宏宇 2022년 9월 8일

편집: 宏宇 2022년 9월 8일

MATLAB Online에서 열기

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 2022년 9월 8일

MATLAB Online에서 열기

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.

댓글을 달려면 로그인하십시오.

Newton iteration doesn't run smoothly

댓글 수: 0
이전 댓글 -2개 표시이전 댓글 -2개 숨기기

답변 (1개)

댓글 수: 2
없음 표시없음 숨기기

참고 항목

카테고리

태그

제품

릴리스

Community Treasure Hunt

Newton iteration doesn't run smoothly

댓글 수: 0 이전 댓글 -2개 표시이전 댓글 -2개 숨기기

답변 (1개)

댓글 수: 2 없음 표시없음 숨기기

참고 항목

카테고리

태그

제품

릴리스

Community Treasure Hunt

댓글 수: 0
이전 댓글 -2개 표시이전 댓글 -2개 숨기기

댓글 수: 2
없음 표시없음 숨기기