How to solve for multiple variables using multiple equations?

I am trying to find the positive values for N1 N2 N3 N4 N5 N6 using the following equations. Running solver gives values for all as 0 by 1. I was suggested to use Newton's iteration but I have no Idea how to do that as I am new to the software.

syms N1 N2 N3 N4 N5 N6;
N7 = 3.135E+00*N1*N2;
N8 = 4.100E-01*N2*N6;
N9 = 1.381E+01*N1*N4;
N10 = 7.610E+01*N2^2*N3;
N11 = 1.720E+05*N2*N3^2*N5^3;
N12 = 4.034E+04*N3*N5^2;
N13 = 2.615E+04*N2*N3^2*N5^2;
N14 = 3.826E+44*N1*N3*N5^2;
N15 = 1.458E+03*N2*N3*N5;
N16 = 2.435E+04*N1*N3^2*N5;
N17 = 1.186E+08*N3^2*N5^3;
N18 = 7.511E+02*N3*N5;
N19 = 6.874E+02*N3*N5;
eqn1 = N1+N2+N3+N4+N5+N6+N7+N8+N9+N10+N11+N12+N13+N14+N15+N16+N17+N18+N19 == 1;
eqn2 = 0.1962*(N2+N7+N8+2*N10+N11+N13+N15)-0.6575*(N1+N7+N9+N14+N16) == 0;
eqn3 = 0.0557*(N3+N10+2*N11+N12+2*N13+N14+N15+2*N16+2*N17+N18+N19)-0.4698*(N4+N9) == 0;
eqn4 = 0.0099*(N5+3*N11+2*N12+2*N13+2*N14+N15+N16+3*N17+N18+N19)-0.3525*(N6+N8) == 0;
eqn5 = 0.0557*(N1+N7+N9+N14+N16)-0.1962*(N4+N9) == 0;
eqn6 = 0.0099*(N2+N7+N8+2*N10+N11+N13+N15)-0.6575*(N6+N8) == 0;
eqns = [eqn1 eqn2 eqn3 eqn4 eqn5 eqn6];
S = solve(eqns,[N1 N2 N3 N4 N5 N6]);

[Stephan]

Is N14 correct with 3.826E+44*N1*N3*N5^2 ? This is much greater then all the other constant factors

[Alex Sha]

if "N14 = 3.826E+44*N1*N3*N5^2 then:

n1: 0.0102966513751017
n2: 0.239882301702705
n3: 0.0756391318035683
n4: 0.0675904371464799
n5: 9.05432958842336E-22
n6: 0.0124929422701635

Feval:
-2.66453525910038E-15
5.39845945723982E-14
-1.38083988687754E-14
1.96308160699266E-12
2.21940521516473E-14
-1.04990321991849E-12

While if "N14 = 3.826E+04*N1*N3*N5^2;", the results will be:

n1: 0.08914189187665
n2: 0.276778705046795
n3: 0.0249154871811356
n4: 0.0266564101251453
n5: 0.00558667307837047
n6: 0.00949311488376945

Feval:
6.66133814775094E-16
1.66533453693773E-16
2.08166817117217E-17
-1.30104260698261E-18
5.20417042793042E-18
7.80625564189563E-18

[TAHMEED TASNIM]

Hello Mr. Alex Sha

Could you please tell me how you were able to get these results? I have the Optimisation toolbox Mr. Alan Stevens suggested, but I dont know how to use it. Is that how you were able to do it?

[Alex Sha]

Hi, the results obtained above are from another math software named 1stOpt, the code looks like below. Matlab should be able to solve this problem, however, it is more possible that some difficulties will be encountered, for example, guess of initial start values. In 1stOpt, guess of start values are no longer needed anymore.

ParameterDomain = [0,];
ConstStr
N7 = 3.135E+00*N1*N2,
N8 = 4.100E-01*N2*N6,
N9 = 1.381E+01*N1*N4,
N10 = 7.610E+01*N2^2*N3,
N11 = 1.720E+05*N2*N3^2*N5^3,
N12 = 4.034E+04*N3*N5^2,
N13 = 2.615E+04*N2*N3^2*N5^2,
N14 = 3.826E+44*N1*N3*N5^2,  //N14 = 3.826E+04*N1*N3*N5^2,
N15 = 1.458E+03*N2*N3*N5,
N16 = 2.435E+04*N1*N3^2*N5,
N17 = 1.186E+08*N3^2*N5^3,
N18 = 7.511E+02*N3*N5,
N19 = 6.874E+02*N3*N5;
Function
N1+N2+N3+N4+N5+N6+N7+N8+N9+N10+N11+N12+N13+N14+N15+N16+N17+N18+N19 = 1;
0.1962*(N2+N7+N8+2*N10+N11+N13+N15)-0.6575*(N1+N7+N9+N14+N16) = 0;
0.0557*(N3+N10+2*N11+N12+2*N13+N14+N15+2*N16+2*N17+N18+N19)-0.4698*(N4+N9) = 0;
0.0099*(N5+3*N11+2*N12+2*N13+2*N14+N15+N16+3*N17+N18+N19)-0.3525*(N6+N8) = 0;
0.0557*(N1+N7+N9+N14+N16)-0.1962*(N4+N9) = 0;
0.0099*(N2+N7+N8+2*N10+N11+N13+N15)-0.6575*(N6+N8) = 0;

Results:

n1: 0.010296651426512
n2: 0.23988230158704
n3: 0.0756391318899229
n4: 0.0675904371072629
n5: 9.05432955846619E-22
n6: 0.0124929422697438

feval:
2.22044604925031E-16
1.0547118733939E-15
2.77555756156289E-17
2.34187669256869E-17
-6.24500451351651E-17
6.24500451351651E-17

Answer 1

[Alan Stevens]

Here's a Newton-Raphson routine for the problem:

% Newton-Raphson Simultaneous Equations
% Six variables N1 to N6 need to be found from the six constraints f1 to f6.
% Variables N7 to N19 are found from N1 to N6.
% 
% The aim is to find values of the N that make each of the f's equal to
% zero (within some tolerance). 

% Set tolerances etc.
f = ones(6,1);  % large initial constraint values
tol = 10^-5;    % constraint tolerance
its = 0;        % iteration counter
flag = true;    % error testing flag

% Initial guesses for N1 to N6
N = [-0.5; 0.5; 0.005; 0.01; 0.2; 0.02];

% Newton-Raphson iteration loop
while flag && its < 60
    
    its = its+1;
    fold = f;
    Nold = N;

    N1 = N(1); N2 = N(2); N3 = N(3); N4 = N(4); N5 = N(5); N6 = N(6);

    k7 = 3.135;
    k8 = 0.41;
    k9 = 13.81;
    k10 = 76.1;
    k11 = 1.720E+05;
    k12 = 4.034E+04;
    k13 = 2.615E+04;
    k14 = 3.826E+04;
    k15 = 1.458E+03;
    k16 = 2.435E+04;
    k17 = 1.186E+08;
    k18 = 7.511E+02;
    k19 = 6.874E+02;

    N7 = k7*N1*N2;            dN7dN1 = k7*N2;         dN7dN2 = k7*N1; 
    N8 = k8*N2*N6;            dN8dN2 = k8*N6;         dN8dN6 = k8*N2; 
    N9 = k9*N1*N4;            dN9dN1 = k9*N4;         dN9dN4 = k9*N1;
    N10 = k10*N2^2*N3;        dN10dN2 = 2*k10*N2*N3;  dN10dN3 = k10*N2^2;
    N11 = k11*N2*N3^2*N5^3;   dN11dN2 = k11*N3^2*N5^3; dN11dN3 = 2*k11*N2*N3*N5^3; dN11dN5 = 3*k11*N2*N3^2*N5^2;
    N12 = k12*N3*N5^2;        dN12dN3 = k12*N5^2;     dN12dN5 = 2*k12*N3*N5;  
    N13 = k13*N2*N3^2*N5^2;   dN13dN2 = k13*N3^2*N5^2; dN13dN3 = 2*k13*N2*N3*N5^2; dN13dN5 = 2*k13*N2*N3^2*N5;
    N14 = k14*N1*N3*N5^2;     dN14dN1 = k14*N3*N5^2;  dN14dN3 = k14*N1*N5^2; dN14dN5 = 2*k14*N1*N3*N5;
    N15 = k15*N2*N3*N5;       dN15dN2 = k15*N3*N5;    dN15dN3 = k15*N2*N5;  dN15dN5 = k15*N2*N3;
    N16 = k16*N1*N3^2*N5;     dN16dN1 = k16*N3^2*N5;  dN16dN3 = 2*k16*N1*N3*N5;   dN16dN5 = k16*N1*N3^2;
    N17 = k17*N3^2*N5^3;      dN17dN3 = 2*k17*N3*N5^3; dN17dN5 = 3*k17*N3^2*N5^2;
    N18 = k18*N3*N5;          dN18dN3 = k18*N5;       dN18dN5 = k18*N3;
    N19 = k19*N3*N5;          dN19dN3 = k19*N5;       dN19dN5 = k19*N3;

    % Jacobian values
    % First row
    J11 = 1 + dN7dN1 + dN9dN1 + dN14dN1 + dN16dN1;
    J12 = 1 + dN7dN2 + dN8dN2 + dN10dN2 + dN11dN2 +dN13dN2 + dN15dN2;
    J13 = 1 + dN10dN3 + dN11dN3 + dN12dN3 + dN13dN3 + dN14dN3 + dN15dN3 + dN16dN3 + dN17dN3 + dN18dN3 + dN19dN3;
    J14 = 1 + dN9dN4;
    J15 = 1 + dN11dN5 + dN12dN5 + dN13dN5 + dN14dN5 + dN15dN5 + dN16dN5 + dN17dN5 + dN18dN5 + dN19dN5;
    J16 = 1 + dN8dN6;
    % Second row
    J21 = 0.1962*dN7dN1 -0.6575*(1 + dN7dN1 + dN9dN1 + dN14dN1 + dN16dN1);
    J22 = 0.1962*(1 + dN7dN2 + dN8dN2 + 2*dN10dN2 + dN11dN2 + dN13dN2 + dN15dN2) -0.6575*dN7dN2;
    J23 = 0.1962*(2*dN10dN3 + dN11dN3 +dN13dN3 + dN15dN3) - 0.6575*(dN14dN3 + dN16dN3);
    J24 = -0.6575*dN9dN4;
    J25 = 0.1962*(dN11dN5 + dN13dN5 + dN15dN5) - 0.6575*(dN14dN5 + dN16dN5);
    J26 = 0.1962*dN8dN6; 
    % Third row
    J31 = 0.0557*(dN14dN1 + 2*dN16dN1) - 0.4698*dN9dN1;
    J32 = 0.0557*(dN10dN2 + 2*dN11dN2 + 2*dN13dN2 + dN15dN2);
    J33 = 0.0557*(1 + dN10dN3 + 2*dN11dN3 + dN12dN3 + 2*dN13dN3 + dN14dN3 +dN15dN3 +2*dN16dN3 + 2*dN17dN3 + dN18dN3 + dN19dN3);
    J34 = -0.4698*(1+dN9dN4);
    J35 = 0.0557*(2*dN11dN5 + dN12dN5 + 2*dN13dN5 + dN14dN5 + dN15dN5 + 2*dN16dN5 + 2*dN17dN5 + dN18dN5 + dN19dN5);
    J36 = 0;
    % Fourth row
    J41 = 0.0099*(2*dN14dN1 + dN16dN1);
    J42 = 0.0099*(3*dN11dN2 + 2*dN13dN2 + dN15dN2) - 0.3525*dN8dN2;
    J43 = 0.0099*(3*dN11dN3 + 2*dN12dN3 +2*dN13dN3 + 2*dN14dN3 + dN15dN3 + dN16dN3 +3*dN17dN3 + dN18dN3 + dN19dN3);
    J44 = 0;
    J45 = 0.0099*(1 + 3*dN11dN5 + 2*dN12dN5 +2*dN14dN5 + dN16dN5 + 3*dN16dN5 + dN18dN5 + dN19dN5);
    J46 = -0.3525*(1 + dN8dN6);
    % Fifth row
    J51 = 0.0557*(1 + dN7dN1 + dN9dN1 + dN14dN1 + dN16dN1) - 0.1962*dN9dN1;
    J52 = 0.0557*dN7dN2;
    J53 = 0.0557*(dN14dN3 + dN16dN3);
    J54 = - 0.1962*(1 + dN9dN4);
    J55 = 0.0557*(dN14dN5 + dN16dN5);
    J56 = 0;
    % Sixth row
    J61 = 0.0099*dN7dN1;
    J62 = 0.0099*(1 + dN7dN2 + dN8dN2 + 2*dN10dN2 + dN11dN2 + dN13dN2 + dN15dN2) - 0.6575*dN8dN2;
    J63 = 0.0099*(2*dN10dN3 + dN11dN3 + dN13dN3 + dN15dN3);
    J64 = 0;
    J65 = 0.0099*(dN11dN5 + dN13dN5 + dN15dN5);
    J66 = 0.0099*dN8dN6 - 0.6575*(1 + dN8dN6);

    % Jacobian matrix
    J = [J11 J12 J13 J14 J15 J16;
         J21 J22 J23 J24 J25 J26;
         J31 J32 J33 J34 J35 J36;
         J41 J42 J43 J44 J45 J46;
         J51 J52 J53 J54 J55 J56;
         J61 J62 J63 J64 J65 J66];

    % Function values
    f1 = N1+N2+N3+N4+N5+N6+N7+N8+N9+N10+N11+N12+N13+N14+N15+N16+N17+N18+N19-1;
    f2 = 0.1962*(N2+N7+N8+2*N10+N11+N13+N15)-0.6575*(N1+N7+N9+N14+N16);
    f3 = .0557*(N3+N10+2*N11+N12+2*N13+N14+N15+2*N16+2*N17+N18+N19)-0.4698*(N4+N9);
    f4 = 0.0099*(N5+3*N11+2*N12+2*N13+2*N14+N15+N16+3*N17+N18+N19)-0.3525*(N6+N8);
    f5 = 0.0557*(N1+N7+N9+N14+N16)-0.1962*(N4+N9);
    f6 = 0.0099*(N2+N7+N8+2*N10+N11+N13+N15)-0.6575*(N6+N8);

    % Function vector
    f = [f1; f2; f3; f4; f5; f6];

    N = Nold - J\f;

    errmax = max(abs(f-fold));
    if errmax>tol
        flag = true;
    else
        flag = false;
    end


end

% Display resulting values of N1 to N6 and constraint values.
% N7 to N19 not displayed here, but appear in the workspace.

str = '       N        f';
disp(str)
disp([N f])

You should check it carefully!

[Alan Stevens]

I should note that in the coefficient of N14 I used 10^04 rather than 10^44.

[TAHMEED TASNIM]

I honestly can't believe that you went through the effort of writing this entire script for a complete stranger.

My hats off to you sir, and thank you very much!!

I am truely grateful for your help.

[Alan Stevens]

I must admit the end result was more extensive than I'd anticipated. However, by the time I realised it would be I was over half way through, so decided to finish rather than waste what I'd done already!

[TAHMEED TASNIM]

I am sorry to ask this after all your help, could you tell me how to make sure all N values are within 0 and 1 in the final result?

I have tried abs and lim([]) functions, but I must be doing something wrong.

[Alan Stevens]

I can't see a legitimate way to do this I'm afraid. It's possible I've not specified all the values of N and their derivatives correctly - you need to check them carefully. An alternative is to use MATLAB's fsolve if you have the Optimisation toolbox (I haven't).