How do I use fsolve with an array input?

I have modelled the three equations below on MATLAB using the fsolve function where N is equal to two. The code associated with the problem is also below. The equations analyse the properties of two components, namely MB and C. There are two types of inputs, scalar and arrays and there are two issues I am facing.

The first issue pertains to the number of solutions I receive. I am expecting to receive four outcomes which would be represented by four columns in the "Solution", however, I only receive three. I initially assumed that because I had three equations, fsolve would only provide me with three solutions, however, upon removing one of the equations, this made no difference to the final result. The three solutions I receive are satisfactory but I just require the fourth (column) now! There are further comments in the code to clarify.

The second issue (less important) relates to writing the sum function. I am attempting to use Fsolve for a summation equation, however, I have no clue with regards to how to formulate this. At the moment, I have simply resorted to using an addition because I only have two components. How would I go about doing this.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
k_MB = 22.213; n_MB = 0.198;
k_C = 35.975;  n_C = 0.234;
%%%%%%%%%%%%%%%%%%%%%%% ARRAY INPUTS %%%%%%%%%%%%%%%%%%%%%%%%%%%%
AC = [3.983 6.558 7.683 5.817]; V = [0.029 0.030 0.032 0.030]; 
Ce_MB = [0.330 0.230 0.188 0.262]; Ce_C = [0.272 0.180 0.141 0.207];
qT = [1.386 0.919 0.868 1.024];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
options=optimset('MaxFunEvals',100000e+30,'MaxIter',100000e+30,'TolFun',.00001);
F = @(x) [(1 - ((((x(3)./((((x(1).*n_MB)./k_MB).^(1/n_MB)) + (x(2).*AC)))))+((x(4)./((((x(1).*n_C)./k_C).^(1/n_C)) + (x(2).*AC))))));...
    (x(1)./x(2)) - (((((x(3)./(((((x(1).*n_MB)./k_MB).^(1/n_MB)) + (x(2).*AC)))*(1/n_C))))+((x(4)./(((((x(1).*n_C)./k_C).^(1/n_C)) + (x(2).*AC))).*(1/n_C)))));...
    ((((x(2) - qT)./qT)+((x(3)-Ce_MB)./Ce_MB)).^2)+ ((((x(2) - qT)./qT)+((x(4)-Ce_C)./Ce_C)).^2)];

% from Function F I hope to derive x(1), x(2), x(3) and x(4), however, at the moment it only provides x(1), x(2) and x(4).

%%%%%%%%%%%%%% SOLUTIONS %%%%%%%%%%%%%%%%%%%%
X0 = [Ce_MB;qT;Ce_C]'; % Initial guess, it won't let me do four initial guesses. 
Solution = fsolve(F,X0)

Spreading_Pressure = Solution(:,1);
qT_Calculated = Solution(:,2);
C_MB = Solution(:,3)
C_C = Solution(:,4); % The function will not provide this as the Solution only has three columns. 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[darova]

I don't understand why your equations are so long? Are summing those value manually?

Also it's unclear which variable is which from the picture

Can you explain?

[Mohammed Rahman]

Thanks for the reply! Yes I am summing them manually which is the issue. I am not sure how to sum them within the fsolve function in a shorter way. I have modelled the variables I am trying to calculate as "x".

Coi = x(3) and x(4)

Phi = x(1)

qT = x(2)

The rest are constants which are defined before the large messy equation. Please let me know if further clarification is needed.

[Torsten]

You need four equations to determine four unknowns, not only three.

[darova]

I don't see this variable on the picture at all: qT = x(2)

Do i have something wrong with my eyes?

[Torsten]

qRj instead of qT

[Torsten]

You obviously set c0i = cRj in your code and remove the absolute values in the expression for S_3 .

Is it this what you want ?

[Mohammed Rahman]

Yes qRJ is x(2) etc. So, my thought was to set the initial guess for cRJ as COI, was there a silly idea? Furthermore, I only have one "S" which I am trying to calculate. I don't quite understand what S_3 means, apologies!

Also your comment about four equations makes sense, however, when I remove one of the equations (leaving two equations and four unknowns) it still outputs three of the four unknowns. Is this a good indication that there is something more fundamental wrong with the code?

Thanks!

[darova]

• Is this a good indication that there is something more fundamental wrong with the code?

Definitely YES. Look at your equations. Are you sure they are correct?

Here is how i'd write your equations:

ftmp1 = phini./Ki).^(1/ni) + qRj*mj./Lj;
f1 = 1 - sum(c0i./ftmp1);
f2 = phi/qRj - sum(c0i./ftmp1./ni);
ftmp2 = abs((qRj-qMj)./qMj) + abs((cRj-cMj)./cMj)
f3 = sum( ftmp.^2 )
F = [f1 f2 f3];

They should be maximum simple and clear!

[Mohammed Rahman]

Thank you much! This solve a big part of our problem! However, I apologise if this is an extremely basic question, but how do I assign the inputs? For example, there are two sets of data (two arrays 2x4) for Coi and two sets of data for ki and ni (two values for k and two values for n). I really do appreciate your help!

C0_MB = [0.330 0.230 0.188 0.262]; 
C0_C = [0.272 0.180 0.141 0.207];

k_MB = 22.213; n_MB = 0.198;
k_C = 35.975;  n_C = 0.234;

[darova]

Ok. i think we are closer to success now

Ki = k_MB and k_C
ni = n_MB and n_C
c0i = C0_MB and C0_C

But who a these? It remains unclear for me

phini = ?
qRj   = ?
qMj   = ?
mj    = ?
Lj    = ?
phi   = ?
cRj   = ?
cMj   = ?

[Mohammed Rahman]

So phini is two separate variables. Phi is one and ni is another. So in the equation it would be phi*ni.

qMj - An array of constant values which are known.

cMj - An array of constant values which are known

mj and Lj = Two constant values which are known.

qRj - One of the unknowns I would like the equations to calculate.

CRj - Two unknown arrays I would like the equations to calculate (for MB and C).

Phi - First term in the denominator of G1. This is an unknown I would like the equations to calculate.

Thanks!

[darova]

Can you fill these constants (put there values)?

phini = ?
qRj   = ?
qMj   = ?
mj    = ?
Lj    = ?
phi   = ?
cRj   = ?
cMj   = ?

[Mohammed Rahman]

%phini = phi*ni
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Ce_MB = [0.330 0.230 0.188 0.262]; 
Ce_C = [0.272 0.180 0.141 0.207];
k_MB = 22.213; 
n_MB = 0.198;
k_C = 35.975;  
n_C = 0.234;
%%%%%%%%% KNOWN %%%%%%%%%%%%%%%%%%%
qMj = [1.386 0.919 0.868 1.024];
mj = [3.983 6.558 7.683 5.817];
Lj = 0.02;
cMj = Ce_MB and Ce_C
ki = k_MB and k_C
ni = n_MB and n_C
%%%%%%%% UNKNOWN %%%%%%%%%%%%%%%%%%
cRj = CR_MB and CR_C % Just to iterate this is two separate arrays.
phi
qRj

[darova]

What is ?

cMj = Ce_MB and Ce_C    % KNOWN
cRj = CR_MB and CR_C    % UKNOWN

[Mohammed Rahman]

Coi is the same as CRj.

Answer 1

[darova]

Here is my attempt without success. No solution found. Is there a mistake i made?

function main
opt = optimset('display','on');
[res_MB,fMB] = fsolve(@(x)func(x,1),ones(1,3),opt);
[res_C,fC] = fsolve(@(x)func(x,2),ones(1,3),opt);
[res_MB(:) fMB(:)]

function res = func(x,ind)
Ce_MB = [0.330 0.230 0.188 0.262]; 
Ce_C = [0.272 0.180 0.141 0.207];
k_MB = 22.213; 
k_C = 35.975;
n_MB = 0.198;
n_C = 0.234;
if ind == 1
    cMj = Ce_MB;
    Ki = k_MB;
    ni = n_MB;
else
    cMj = Ce_C;
    Ki = k_C;
    ni = n_C;
end
qMj = [1.386 0.919 0.868 1.024];
mj = [3.983 6.558 7.683 5.817];
Lj = 0.02;

%%%%%%%% UNKNOWN %%%%%%%%%%%%%%%%%%
% cRj = CR_MB and CR_C % Just to iterate this is two separate arrays.
% phi
% qRj

phi = x(1);
qRj = x(2);
cRj = x(3);
c0i = x(3);
ftmp1 = (phi.*ni./Ki).^(1/ni) + qRj*mj./Lj;
f1 = 1 - sum(c0i./ftmp1);
f2 = phi/qRj - sum(c0i./ftmp1./ni);
ftmp2 = abs((qRj-qMj)./qMj) + abs((cRj-cMj)./cMj);
f3 = sum( ftmp2.^2 );
res = [f1 f2 f3];

[Torsten]

f3 can never be 0 if there are different values for the qMj or the cMj.

And if they are equal, qRj must equal qMj and cRj must equal cMj.

[Mohammed Rahman]

Hi Torsten, thanks for the reply! So the point in F3 is just to evaluate the difference between the calculated values (from f1 and f2) and the pre-defined values (cMj and qMJ). I am not expecting F3 to be equal to zero. My Classmate and I initially attempted to solve the system using Newton Raphson to no avail which involved minising f3.

Thank you for the answer Darova!! My classmate and I are going through your code, I don't know how to repay you!

[darova]

• I don't know how to repay you!

But i do!

[darova]

Here is another way to analyze your equations: visually

x1 = linspace(-10,10,20);
y1 = linspace(-1,1,20);
% phi = x
% qRj = y
% cRj = z
% c0i = z
[x,y,z] = meshgrid(x1,y1,x1);
G1 = x*0;
G2 = x*0;
S3 = x*0;
for i = 1:size(x,1)
    for j = 1:size(x,2)
        for k = 1:size(x,3)
            xx = [x(i,j,k) y(i,j,k) z(i,j,k)];
            res = func(xx,2);
            G1(i,j,k) = res(1);
            G2(i,j,k) = res(2);
            S3(i,j,k) = res(3);
        end
    end
end
cla
p1 = isosurface(x,y,z,G1,0);
p2 = isosurface(x,y,z,G2,0);
p3 = isosurface(x,y,z,S3,20);
patch(p1,'facecolor','r')
patch(p2,'facecolor','g')
patch(p3,'facecolor','b')
light
axis vis3d 
% alpha(0.5)
legend('G1==0','G2==0','S3==20')

[Mohammed Rahman]

Hi Darova,

Apologies for the late follow up, however, my classmate and I were wondering what the output (ans) corresponds to? I.e which variables does the code output in the answer you sent?

Thanks!

[darova]

You mean this?

% res_MB are values phi qRj cRj/c0i
% fMb are results function: G1 G2 S3
[res_MB(:) fMB(:)]

See more fsolve

[Mohammed Rahman]

Yes! Thank you so much, that makes sense now.

[darova]

you are welcome

[Mohammed Rahman]

Hi Darova!

Hope you're doing well in this state of quarantine. I just had a quick follow up, my classmate and I have been slaving away at this code and have simplified it further but have run into a tiny issue. We have boiled down the problem and no we simply have two equations with two unknowns; phi and qRj. Our issue is, when we run the code we obtain two 2x2 matrices as answers whereas we would like one 2x2 matrix. We believe the code is solving f1 and f2 for the two sets of data (DATA SET 1 and DATA SET 2) separately rather than actually summing it.

How would we adapt this code to return one 2x2 matrix which contains a value for phi, qRj, f1 and f2?

THANKS!

function IAST
opt = optimset('display','on');
x0 = ones(1,2);
[res_1,f1] = fsolve(@(x)func(x,1),x0,opt);
[res_2,fc2]= fsolve(@(x)func(x,2),x0,opt);
[res_1(:) f1(:)]
[res_2(:) fc2(:)]

function res = func(x,ind)
%%%%%%%% DATA SET 1
C0i_1 = 5.651086447;
k_1 = 35.97493352;
n_1 = 0.2338;

%%%%%%% DATA SET 2
C0i_2 = 9.488130182;
k_2 = 22.21263089;
n_2 = 0.1981;

%%%%%%%%%%%%%%%%%%%%
mj = 11.506;

Lj = 0.02;

if ind == 1
    c0i = C0i_1;     %%% Initial Concentration - Component 1
    Ki = k_1;        %%% Freundlich K - Component 1
    ni = n_1;        %%% Freundlich n - Component 1
else
    c0i = C0i_2;     %%% Initial Concentration - Component 2
    Ki = k_2;        %%% Freundlich K - Component 2
    ni = n_2;        %%% Freundlich n - Component 2
end

phi = x(1);
qRj = x(2);
%%%%%%%%%%%%%%%%%%%%%%%%% EQUATIONS %%%%%%%%%%%%%%%%%%%%%%%%%%
ftmp1 = (phi*ni/Ki)^(1/ni) + qRj*(mj/Lj);
f1 = 1 - sum(c0i/ftmp1);
f2 = phi/qRj - sum(c0i/ftmp1/ni);
res = [f1 f2];

[darova]

res_1 - roots after first solving equations, f1 values of function at those points (roots)

res_2 - roots after second solving equations, fc2 values of function at those points (roots)

[res_1(:) f1(:)]
[res_2(:) fc2(:)]

So 8 values together

[Mohammed Rahman]

Ahhh this makes a lot more sense now! My original intention was to have a value of phi and crj which satisfy both equations; rather than having the roots of both functions (essentially a simultaneous equation problem). Is fsolve the wrong function to carry this out?

Just an update - I attempted to add the func(x,2) to the resolution function, does this solve my problem above or is this completely wrong? Btw I really appreciate the help, cannot thank you enough.

function IAST
opt = optimset('display','on');
x0 = ones(1,2);
[res_1,f1] = fsolve(@(x)func(x,1)+func(x,2),x0,opt); %%% Does this make sense?
[res_2,fc2]= fsolve(@(x)func(x,2),x0,opt);
[res_1(:) f1(:)]
[res_2(:) fc2(:)]

[darova]

It's impossible, for each set of data you will have unique solution

c1x + c2y = 0
c3x - c3y = 0

SOlution: x1 and y1

c5x + c6y = 0
c7x - c8y = 0

SOlution: x2 and y2

But you can't make x1=x2 and y1=y2

[Mohammed Rahman]

Thanks for the reply! I see what you're saying. Our system sums the two sets of data in f1 and f2, so in essence I have understood the system to be what is shown below, so in theory I should two equations in the system. The nifty notation you showed me simplfied the huge equation I showed initially which added together all the terms separately. Or is my understanding fundamentally incorrect?

c1x + c2y + c5x + c6y = 0
c3x - c3y + c7x - c8y = 0

ftmp1 = (phi*ni/Ki)^(1/ni) + qRj*(mj/Lj);
f1 = 1 - sum(c0i/ftmp1); % Sum which contains all the constants from the data sets.
f2 = phi/qRj - sum(c0i/ftmp1/ni); % Sum which contains all the constants from the data sets.

[darova]

[Mohammed Rahman]

There are only two unknowns (x and y), highlighted below in yellow and everything else is a constant (c1..cn). The two large brackets are what the sum equation is trying to represent.

[darova]

Im lost. Dont know what to say. Is this another question or the same?

[Mohammed Rahman]

Apologies for the confusion! We have figured it out, it was the same problem. We really appreciate your help!!

[darova]

you are welcome

[Mohammed Rahman]

Hi Darova,

Apologies about popping up again, however, I was wondering if you could help me! This is another variation of your answer, however, my classmate and I are having massive issues with regards to the variable mj (Line 4). In our solve function; Phi_qRJ (Line 15), we would like to solve for each element of the mj array, is this possible? Whenever we try to use the arrayfun with fsolve we keep getting "matrix dimesions do not not agree".

mj is used in:

• Line 6 (qAC)
• Line 27 (Z2)
• Line 28 (Z3)

In essence I would like the code to use one value of mj at a time, rather than using the entire array. I appeciate any help you can offer!

function IAST_2
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% INPUTS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
global mj Lj qAC K n C0i Phi qRj
mj = [0.1,0.2,0.3,0.4,0.5];      % Mass of Activated Carbon (g)
Lj = 0.1;                        % Volume of Solution (dm3)
qAC = mj./Lj;                    % Activated Carbon Loading in Solution (g/dm3)

K = [49.1103, 69.49];
n = [0.2095, 0.19727];
C0i = [1,1];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% SOLVE %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
opt = optimset('display','on'); 
x0 = ones(1,2);

Phi_qRj = fsolve(@(x)G1_G2(x),x0,opt);

Phi = Phi_qRj(:,1)
qRj = Phi_qRj(:,2)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function Equations_G1_G2 = G1_G2(x)
global mj Lj qAC K n C0i Phi qRj
Phi = x(1);
qRj = x(2);

Z1 = ((Phi*n)./K).^(1./n)
Z2 = (qRj * qAC)
Z3 = C0i./(Z1+Z2);

Z4 = C0i./(Z1+Z2).*(1./n);

G1 = sum(Z3) + -1;
G2 = sum(Z4) + -(Phi/qRj);

Equations_G1_G2 = [G1 G2]';

[darova]

Try for loop