Newton Raphson Liquid-liquid extraction

Hello everyone, I wonder if someone could help me resolving an error Im having on my matlab code. I am making a program based on a liquid liquid extraction. I made a function with 8 equations and 8 variables for finding de molar fraction x. The error matlab is reporting correspond to “index exceeds the number of array elements” I don’t know what I have left to do. Here is my code

clc
clear
close all
n = 8;
x1space = linspace(0,1,n);
xresp = zeros(1,n);
Tresp = zeros(1,n);
for i = 1:n
    iter = 0;
    iterMax = 1000;
    error = 1;
    tol = 1e-6;
    
    x0 = (0:n:1)
    while error > tol && iter < iterMax
        fx0 = BM(x0);
        jx0 = jacobian(x0);
        delta = -jx0\fx0;
        xnew = x0 + delta;
        error = norm(delta);
        x0 = xnew;
    end
    xresp(i) = xnew(1:8);
    
end

x0 = 0

Index exceeds the number of array elements. Index must not exceed 1.

Error in solution>BM (line 30)
x1 = u(1) ; x2 = u(2); x3= u(3); x4= u(4); x5 = u(5) ; x6 = u(6); x7 = u(7); x8 = u(8);

function sol = BM(u)
WA = 100; 
WS = 10; 
x0 = 0.05;
x1 = u(1) ; x2 = u(2); x3= u(3); x4= u(4); x5 = u(5) ; x6 = u(6); x7 = u(7); x8 = u(8); 
sol1= WA*x0 +WS*(12*x2+(18*x2)^2)-(WA*x1+WS*(12*x1+18*(x1)^2)); 
sol2= WA*x1+WS*(12*x3+(18*x3)^2)-(WA*x2+WS*(12*x2+18*(x2)^2)); 
sol3= WA*x2+WS*(12*x4+(18*x4)^2)-(WA*x3+WS*(12*x3+18*(x3)^2)); 
sol4= WA*x3+WS*(12*x5+(18*x5)^2)-(WA*x4+WS*(12*x4+18*(x4)^2)); 
sol5= WA*x4 *WS*(12*x6+(18*x6)^2)-(WA*x5+WS*(12*x5+18*(x5)^2)); 
sol6= WA*x5 *WS*(12*x7+(18*x7)^2)-(WA*x6+WS*(12*x6+18*(x6)^2)); 
sol7= WA*x6 *WS*(12*x8+(18*x8)^2)-(WA*x7+WS*(12*x7+18*(x7)^2)); 
sol8= WA*x7 *WS*(12*x0+(18*x0)^2)-(WA*x8+WS*(12*x8+18*(x8)^2)); 
end 
function j = jacobian(u)
y = BM(u);
h = 1e-8;
for i=1:numel(u)
    u(i) = u(i) + h;
    yh = BM(u);
    j(:,i) = (yh-y)/h;
    u(i) = u(i) - h;
end 
end

I’ll be really grateful if someone could help me to resolve it.

[Walter Roberson]

Does this question mean that you have solved https://www.mathworks.com/matlabcentral/answers/1674334-newton-raphson-multivariate-liquid-liquid-extraction?s_tid=srchtitle ?

Answer 1

[Walter Roberson]

n = 8

n = 8

x0 = (0:n:1)

x0 = 0

That is a vector of length 1, because it is "start at 0, add 8 at a time, stop if you exceed 1".

Your code requires that it is a vector of length 8.

Hint:

0:n-1

Answer 2

[Torsten]

No sum constraint on sum_i xi ?

n = 8;
x1space = linspace(0,1,n);
xresp = zeros(1,n);
Tresp = zeros(1,n);
%for i = 1:n
iter = 0;
iterMax = 1000;
error = 1;
tol = 1e-6;

x0 = 1/n*ones(n,1);
while error > tol && iter < iterMax
    fx0 = BM(x0);
    jx0 = jacobian(x0);
    delta = -jx0\fx0;
    xnew = x0 + delta;
    error = norm(delta);
    x0 = xnew;
end
xresp = xnew
BM(xresp)

%end
end
function sol = BM(u)
WA = 100; 
WS = 10; 
x0 = 0.05;
x1 = u(1) ; x2 = u(2); x3= u(3); x4= u(4); x5 = u(5) ; x6 = u(6); x7 = u(7); x8 = u(8); 
sol = zeros(8,1);
sol(1)= WA*x0 +WS*(12*x2+(18*x2)^2)-(WA*x1+WS*(12*x1+18*(x1)^2)); 
sol(2)= WA*x1+WS*(12*x3+(18*x3)^2)-(WA*x2+WS*(12*x2+18*(x2)^2)); 
sol(3)= WA*x2+WS*(12*x4+(18*x4)^2)-(WA*x3+WS*(12*x3+18*(x3)^2)); 
sol(4)= WA*x3+WS*(12*x5+(18*x5)^2)-(WA*x4+WS*(12*x4+18*(x4)^2)); 
sol(5)= WA*x4 *WS*(12*x6+(18*x6)^2)-(WA*x5+WS*(12*x5+18*(x5)^2)); 
sol(6)= WA*x5 *WS*(12*x7+(18*x7)^2)-(WA*x6+WS*(12*x6+18*(x6)^2)); 
sol(7)= WA*x6 *WS*(12*x8+(18*x8)^2)-(WA*x7+WS*(12*x7+18*(x7)^2)); 
sol(8)= WA*x7 *WS*(12*x0+(18*x0)^2)-(WA*x8+WS*(12*x8+18*(x8)^2)); 
end 
function j = jacobian(u)
y = BM(u);
h = 1e-8;
for i=1:numel(u)
    u(i) = u(i) + h;
    yh = BM(u);
    j(:,i) = (yh-y)/h;
    u(i) = u(i) - h;
end 
end

Answer 3

[Juliana Quintana]

Thank you for the answers. I run the program and unfortunately the answer found doesn’t seems right

[Torsten]

As last time: Maybe you'll find an error in your equations. E.g. I don't see a sum constraint on the x_i, meaning that you force them to sum to 1 and in sol(5) to sol(8), you write

sol(5)= WA*x4 *WS*(12*x6+(18*x6)^2)-(WA*x5+WS*(12*x5+18*(x5)^2)); 
sol(6)= WA*x5 *WS*(12*x7+(18*x7)^2)-(WA*x6+WS*(12*x6+18*(x6)^2)); 
sol(7)= WA*x6 *WS*(12*x8+(18*x8)^2)-(WA*x7+WS*(12*x7+18*(x7)^2)); 
sol(8)= WA*x7 *WS*(12*x0+(18*x0)^2)-(WA*x8+WS*(12*x8+18*(x8)^2)); 

instead of

sol(5)= WA*x4 +WS*(12*x6+(18*x6)^2)-(WA*x5+WS*(12*x5+18*(x5)^2)); 
sol(6)= WA*x5 +WS*(12*x7+(18*x7)^2)-(WA*x6+WS*(12*x6+18*(x6)^2)); 
sol(7)= WA*x6 +WS*(12*x8+(18*x8)^2)-(WA*x7+WS*(12*x7+18*(x7)^2)); 
sol(8)= WA*x7 +WS*(12*x0+(18*x0)^2)-(WA*x8+WS*(12*x8+18*(x8)^2));