MATLAB: How to solve for multiple variables using multiple equations

newton iterationsolver

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]);

Best Answer

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!
Related Question