How do I fix this error?

Hi, I want to solve a set of equations of 15 equations and 15 unknowns using Newton-Raphson method, but after writing the program and executing it gives the following error:

Error using symengine

Arithmetical expression expected.

Error in sym/privUnaryOp (line 1019)

Csym = mupadmex(op,args{1}.s,varargin{:});

Error in sym/abs (line 9)

Y = privUnaryOp(X, 'symobj::map', 'abs');

Error in new (line 23)

if abs(xold-xnew)<1e-14

How do I fix these errors?

Program text:

My function by name newfunc:

function EQN=newfun(x,P,R,phi,xold)

Y_CO2=x(1)    ; Y_H2O=x(2)    ; Y_N2=x(3)     ; Y_CO=x(4) ; Y_O2=x(5)  ;
Y_CH3OH=x(6)  ; Y_CH2O=x(7)   ; Y_CH4=x(8)    ; Y_H2=x(9) ;
Landa_C=x(10) ; Landa_H=x(11) ; Landa_O=x(12) ; 
Landa_N=x(13) ; Landa_sigma=x(14) ; T=x(15);

Gs.CO2=-393.360+((-3.8212e-03)*T)+(1.3322e-06)*T^2 ; % from article
Gs.H2O=((2.053e-06)*T^2)+(0.0496*T)-2.436e+02 ;  %unit: J/Mol
Gs.N2=0 ;
Gs.CO=-109.885+(-9.2218E-02*T)+1.4547E-06*T^2 ;% from article
Gs.O2=0 ;
Gs.CH3OH=-201.860+(1.2542e-01*T)+(2.0345e-05)*T^2 ;% from article
Gs.CH2O=((1.3642e-06)*T^2)+(0.0333*T)-1.216658e+02 ;%from curve fitting
Gs.CH4=-201.860 +(1.2542e-01*T)+((2.0345e-05)*T^2) ;% from article
Gs.H2=0 ;

Cp_CO2=27.437+((4.2315e-02)*T)+((-1.9555e-05)*T^2)+((+3.9968e-09)*T^3)+(-2.9872e-13)*T^4;
Cp_H2O=33.933+((-8.4186e-03)*T)+((+2.9906e-05)*T^2)+((-1.7825e-08)*T^3)+(3.6934e-12)*T^4;
Cp_N2=29.342+((-3.5395e-03)*T)+((+1.0076e-05)*T^2)+((-4.3116e-09)*T^3)+(2.5935e-13)*T^4;
Cp_CO=29.556+((-6.5807e-03)*T)+(+2.0130e-05*T^2)+(-1.2227e-08*T^3)+(2.2617e-12)*T^4;
Cp_O2=29.526+((-8.8999e-03)*T)+((+3.8083e-05)*T^2)+((-3.2629e-08)*T^3)+(8.8607e-12)*T^4;
Cp_CH3OH=40.046+((-3.8287e-02)*T)+((2.4529e-04)*T^2)+((-2.1679e-07)*T^3)+(5.9909e-11)*T^4;
Cp_CH2O=((-5.1535e-13)*T^4)+((6.5842e-09)*T^3)+((-3.1930e-05)*T^2)+(0.0717*T)+15.856;%from curve fitting
Cp_CH4=34.942+((-3.9957e-02)*T)+((1.9184e-04)*T^2)+((-1.5303e-07)*T^3)+(3.9321e-11)*T^4;
Cp_H2=25.399+((2.0178e-02)*T)+((-3.8549e-05)*T^2)+((3.1880e-08)*T^3)+(-8.7585e-12)*T^4;
delta_Hp=Y_CO2*int(Cp_CO2,T,298,T)+ Y_H2O*int(Cp_H2O,T,298,T)+...
    Y_CO*int(Cp_CO,T,298,T)+Y_N2*int(Cp_N2,T,298,T)+...
    Y_O2*int(Cp_O2,T,298,T)+Y_CH3OH*int(Cp_CH3OH,T,298,T)+...
    Y_CH2O*int(Cp_CH2O,T,298,T)+Y_CH4*int(Cp_CH4,T,298,T)+...
    Y_H2*int(Cp_H2,T,298,T);
delta_Hs = Y_CO2*-393509+ Y_H2O*-241818+ Y_CO*-110525+...
    Y_CH3OH*-200660+ Y_CH2O*-108570+ Y_CH4*-74520-(-238660/Landa_sigma); % delta_H_298

eqn(1)=Gs.CO2+R*T*taylor(log(P*phi*Y_CO2),Y_CO2,'ExpansionPoint',abs(xold(1)))+Landa_C+2*Landa_O;
eqn(2)=Gs.H2O+R*T*taylor(log(P*phi*Y_H2O),Y_H2O,'ExpansionPoint',abs(xold(2)))+2*Landa_H+Landa_O;
eqn(3)=Gs.N2+R*T*taylor(log(P*phi*Y_N2),Y_N2,'ExpansionPoint',abs(xold(3)))+2*Landa_N;
eqn(4)=Gs.CO+R*T*taylor(log(P*phi*Y_CO),Y_CO,'ExpansionPoint',abs(xold(4)))+Landa_C+Landa_O;
eqn(5)=Gs.O2+R*T*taylor(log(P*phi*Y_O2),Y_O2,'ExpansionPoint',abs(xold(5)))+2*Landa_O;
eqn(6)=Gs.CH3OH+R*T*taylor(log(P*phi*Y_CH3OH),Y_CH3OH,'ExpansionPoint',abs(xold(6)))+Landa_O+4*Landa_H+Landa_C;
eqn(7)=Gs.CH2O+R*T*taylor(log(P*phi*Y_CH2O),Y_CH2O,'ExpansionPoint',abs(xold(7)))+Landa_O+2*Landa_H+Landa_C;
eqn(8)=Gs.CH4+R*T*taylor(log(P*phi*Y_CH4),Y_CH4,'ExpansionPoint',abs(xold(8)))+4*Landa_H+Landa_C;
eqn(9)=Gs.H2+R*T*taylor(log(P*phi*Y_H2),Y_H2,'ExpansionPoint',abs(xold(9)))+2*Landa_H;
eqn(10)=Y_CO2+Y_CO+Y_CH3OH+Y_CH2O+Y_CH4-(1/Landa_sigma); %k=C
eqn(11)=2*Y_H2O+4*Y_CH3OH+2*Y_CH2O+4*Y_CH4+2*Y_H2-(4/Landa_sigma); %k=H
eqn(12)=2*Y_CO2+Y_H2O+Y_CO+2*Y_O2+Y_CH3OH+Y_CH2O-(4/Landa_sigma); %k=O
eqn(13)=2*Y_N2-(11.28/Landa_sigma); %k=N
eqn(14)=Y_CO2+Y_H2O+Y_N2+Y_CO+Y_O2+Y_CH3OH+Y_CH2O+Y_CH4+Y_H2-1; 
eqn(15)=delta_Hs + delta_Hp; % for calculation Tmax
EQN=eqn;
end

Execution text with name new:

clc;clear all;close all;
%% Inputs
P=1;%input('Pressure(atm)= ');
R=0.008314;%unit: KJ/Mol.K
phi=0.5:0.1:1.5;%input('PHI= ');
n=15;%input('Number of unknown= ');
xold=ones(1,n)*1e-03;%input('Initial guess= ');
disp(' ');
if length(xold)~=n
    error('The dimention of initial guess are not correct');
end
%% Calculation of newton rafson
xold=xold';
x=sym('x',[1 n]);
for j=1:length(phi)
    for i=1:1000
        F=vpa(newfun(x,P,R,phi(j),xold),14);
        J=vpa(jacobian(F,x),14);%vpa:Display the value of the sub with 14 decimal places
        E_F=vpa(subs(F,x,xold'),14);%sub: Replace the parameter
        E_J=vpa(subs(J,x,xold'),14);
        delta_x=vpa(inv(E_J)*-E_F',14);
        xnew=vpa(xold+delta_x,14);
        if abs(xold-xnew)<1e-14
            break
        end
        if i==1000
            error('Maximum iteration')
        end
        xold=xnew;
    end
    j
    Result(:,j)=xnew;
end
y_total=sum(xnew(1:9));
%% Outputs
disp('"Results"') ; disp(' ')     ; disp('Y_Total= ') ; disp(y_total);
disp('Y_CO2= ')   ; disp(xnew(1)) ; disp('Y_H2O= ')   ; disp(xnew(2));
disp('Y_N2= ')    ; disp(xnew(3)) ; disp('Y_CO= ')    ; disp(xnew(4));
disp('Y_O2= ')    ; disp(xnew(5)) ; disp('Y_CH3OH= ')    ; disp(xnew(6));
disp('Y_CH2O= ')  ; disp(xnew(7)) ; disp('Y_CH4= ')     ; disp(xnew(8));
disp('Y_H2= ')    ; disp(xnew(9)) ;
disp('Landa_C= ') ; disp(xnew(10));
disp('Landa_H= ') ; disp(xnew(11)); disp('Landa_O= ') ; disp(xnew(12));
disp('Landa_N= ') ; disp(xnew(13)); disp('Landa_sigma= ');disp(xnew(14));
disp('T=')        ; disp(xnew(15));

[Walter Roberson]

By the time i == 5 for j == 1, then

        E_J=vpa(subs(J,x,xold'),14);

that matrix becomes singular, which means there is no solution for the system using this approach.

• perhaps newfun() has some mistake
• perhaps your Euler implementation is wrong
• perhaps there is no solution for the initial conditions

[fatemeh torbati]

thank you so much...