Symbolic parameters not supported in nonpolynomial equations. in Newton iteration

The purpose of this program is to solve 'x' by Newton iteration.

The author knows that there is a problem with 'vpasolve', but there is still a problem after changing to 'solve'.

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.

% Set T2 to x.
function F=funfqtotceshi7(x)
global R q_tot fq_tot y i A B C TM K E S
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=vpasolve((C1P*a+C2f*(0.017+(7E-6)*(800-(x+xn1)/2)+0.0228*log((x+xn1)/2))/D)*(x-xn1)+(x.^4-xn1.^4)*g/(1/c+1/c-1)-q_tot==0,xn1);
eq2=vpasolve((C1P*a+C2f*(0.017+(7E-6)*(800-(eq1+xn2)/2)+0.0228*log((eq1+xn2)/2))/D)*(eq1-xn2)+(eq1.^4-xn2.^4)*g/(1/c+1/c-1)-q_tot==0,xn2);
eq3=vpasolve((C1P*a+C2f*(0.017+(7E-6)*(800-(eq2+xn3)/2)+0.0228*log((eq2+xn3)/2))/D)*(eq2-xn3)+(eq2.^4-xn3.^4)*g/(1/c+1/c-1)-q_tot==0,xn3);
eq4=vpasolve((C1P*a+C2f*(0.017+(7E-6)*(800-(eq3+xn4)/2)+0.0228*log((eq3+xn4)/2))/D)*(eq3-xn4)+(eq3.^4-xn4.^4)*g/(1/c+1/c-1)-q_tot==0,xn4);
eq5=vpasolve((C1P*a+C2f*(0.017+(7E-6)*(800-(eq4+xn5)/2)+0.0228*log((eq4+xn5)/2))/D)*(eq4-xn5)+(eq4.^4-xn5.^4)*g/(1/c+1/c-1)-q_tot==0,xn5);
eq6=vpasolve((C1P*a+C2f*(0.017+(7E-6)*(800-(eq5+xn6)/2)+0.0228*log((eq5+xn6)/2))/D)*(eq5-xn6)+(eq5.^4-xn6.^4)*g/(1/c+1/c-1)-q_tot==0,xn6);
A=[T1 x eq1 eq2 eq3 eq4 eq5 eq6];
%the first round of temperature
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=1:7
    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 i=1:7
    y(1)=TC;
    y(i+1)=TC+((sum(R(8-i: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

Differentiate Equation dF

function dF=dfunfqtotceshi4(x)
global R q_tot fq_tot y i A B C TM K E S 
syms x
T1=300; % assumption 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'
syms xn1 xn2 xn3 xn4 xn5 xn6
eq1=vpasolve((C1P*a+C2f*(0.017+(7E-6)*(800-(x+xn1)/2)+0.0228*log((x+xn1)/2))/D)*(x-xn1)+(x.^4-xn1.^4)*g/(1/c+1/c-1)-q_tot==0,xn1);
eq2=vpasolve((C1P*a+C2f*(0.017+(7E-6)*(800-(eq1+xn2)/2)+0.0228*log((eq1+xn2)/2))/D)*(eq1-xn2)+(eq1.^4-xn2.^4)*g/(1/c+1/c-1)-q_tot==0,xn2);
eq3=vpasolve((C1P*a+C2f*(0.017+(7E-6)*(800-(eq2+xn3)/2)+0.0228*log((eq2+xn3)/2))/D)*(eq2-xn3)+(eq2.^4-xn3.^4)*g/(1/c+1/c-1)-q_tot==0,xn3);
eq4=vpasolve((C1P*a+C2f*(0.017+(7E-6)*(800-(eq3+xn4)/2)+0.0228*log((eq3+xn4)/2))/D)*(eq3-xn4)+(eq3.^4-xn4.^4)*g/(1/c+1/c-1)-q_tot==0,xn4);
eq5=vpasolve((C1P*a+C2f*(0.017+(7E-6)*(800-(eq4+xn5)/2)+0.0228*log((eq4+xn5)/2))/D)*(eq4-xn5)+(eq4.^4-xn5.^4)*g/(1/c+1/c-1)-q_tot==0,xn5);
eq6=vpasolve((C1P*a+C2f*(0.017+(7E-6)*(800-(eq5+xn6)/2)+0.0228*log((eq5+xn6)/2))/D)*(eq5-xn6)+(eq5.^4-xn6.^4)*g/(1/c+1/c-1)-q_tot==0,xn6);
A=[T1 x eq1 eq2 eq3 eq4 eq5 eq6];
%the first round of temperature
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=1:7
    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 i=1:7
    y(1)=TC;
    y(i+1)=TC+((sum(R(8-i: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); %Take the derivative of x in F.
end

Newton iteration

% Newton iteration
clear;clc
syms x
F=funfqtotceshi7(x);
dF=dfunfqtotceshi4(x);
x(1)=290; % Define the iteration initial value.
for ii=1:10000 % iterations
    dt(ii)=funfqtotceshi7(x(ii))./dfunfqtotceshi4(x(ii));
    x(ii+1)=x(ii)-dt(ii);
    if abs(dt(ii))~=0.001
        break;
    else
    end
    x0=x; % update iterative result
end
disp('Solve results：');
x(ii+1)
disp('number of iterations：');
ii

Many thanks in adavance,

Hopefully there is a solution for this problem.

Answer 1

[Nipun]

Hi,

I understand that you are working on a Newton iteration program to solve for 'x' in a set of equations, and you are facing some issues. Having gone through your code, I have identified potential problems. Consider the following recommendations to help with the errors.

1. Replace vpasolve with solve when solving the system of equations. vpasolve is used for numeric solutions, and it might not be necessary in this case. The link to the documentation for solve function is given at the end of the answer
2. Avoid redefining the symbol x in the dfunfqtotceshi4 function since it's already defined as the argument of the function.
3. Instead of checking for abs(dt(ii)) ~= 0.001, consider using a more meaningful stopping criterion. For example, you could check if abs(dt(ii)) is below a certain tolerance, like 1e-6.
4. Consider preallocating the arrays to improve performance.

I am attaching the mended code based on the above recommendations, below:

function F = funfqtotceshi7(x)
    % ... (unchanged)
end

function dF = dfunfqtotceshi4(x)
    % ... (unchanged)
end

% Newton iteration
clear; clc
syms x
F = funfqtotceshi7(x);
dF = dfunfqtotceshi4(x);

x(1) = 290; % Initial guess
maxIterations = 10000;
tolerance = 1e-6; % Set a reasonable tolerance

dt = zeros(1, maxIterations);

for ii = 1:maxIterations
    dt(ii) = F / dF;
    x(ii + 1) = x(ii) - dt(ii);

    if abs(dt(ii)) < tolerance
        disp('Converged.');
        break;
    end
end

disp('Solve results:');
disp(x(ii + 1));
disp('Number of iterations:');
disp(ii);

Link to documentation:

1. Solve function: https://www.mathworks.com/help/symbolic/sym.solve.html

Hope this helps.

Regards,

Nipun