How to solve for N sets of 2 equations and 2 unknowns using a for loop?

Hello everyone!

I need to solve a problem involving N number of 2-equations and 2-unknowns inside a for-loop. The N sets of equations share simiar structure thus running them in a for loop.

for example (values are not check, they are presented just to show illustration of the problem):

N = 3;

a = [1,2,3];

b = [2,1,3];

d = [3,2,2];

e = [2,2,1];

C1 = [ 10, 4, 6];

C2 = [ 15, 3, 4];

for i = 1:N

a(i)*x^2+b(i)*y = C1(i);

d(i)*x+e(i)*y = C2(i);

fsolve/vpasolve

end

Answer 1

[Jon]

This is one approach using nested functions to pass additional paramaters, see https://www.mathworks.com/help/optim/ug/passing-extra-parameters.html

Note: only the second set of parameters returned an exit flag of 1, solutions were not obtained for the first and last sets. You would have to look into this more deeply if these are really the values you are using. Maybe you need a better initial guess, or maybe there really isn't a solution for that set of values.

Also, in general when solving systems of non-linear equations there may be more than one solution. If this is the case you will need to do more work to make sure you get the one you want, fsolve will just pass back the first solution it finds.

% define vectors with coefficient (parameter) values
a = [1,2,3];
b = [2,1,3];
d = [3,2,2];
e = [2,2,1];
C1 = [ 10, 4, 6];
C2 = [ 15, 3, 4];

% determine number of evaluations needed
N = numel(a); % assume they all have same length

% preallocate arrays to hold results
X = zeros(2,N);
exitflags = zeros(1,N)
% assign initial guess, assume same for all loops
x0 = zeros(2,1);

% loop to solve for each set of parameters
for k = 1:N % use k as loop counter, MATLAB uses i for imaginary values
    % you should check exitflag values, should be 1 for reliable answer
    [X(:,k),exitflags(k)] = solveSys(a(k),b(k),d(k),e(k),C1(k),C2(k),x0);
end

% use nested functions to pass additional parameters
function [X,exitflag] = solveSys(a,b,d,e,C1,C2,x0)
[X,~,exitflag] = fsolve(@fun,x0);

    function fval = fun(x)
    % define function which will equal zero at desired values of x
    % Note: your x = x(1) and your y = x(2)
    % since this is a nested function the values of a, b, ... and
    % loop counter from main program
    % will be available here
    fval(1) = a*x(1)^2 + b*x(2) - C1 ;
    fval(2) = d*x(1) + e*x(2) -C2;
    end
end

[John Harvie delos Reyes]

Thank you very much!

I have a follow up question tho, is it possible to have the values of the initial guess run inside the loop (every initial guess is different per loop)?

Edit: I figured it out. will work on my original problem now. Thank you very much once again!

[John Harvie delos Reyes]

Hello!

I finished writing my code and it runs no error but is not able to provide me meaningful answers. Can you please take a look at my code?

I might have missed something while im translating the example into my actual problem.

clear;

clc;

%%

Disp = [1 10 4 2 3;

2 8 3 2 2

3 5 4 1.5 3];

options = optimset('TolFun',1e-12,'MaxFunEvals',1e9,'Maxiter',1e9);

options.StepTolerance =1e-9;

%%

N = size(Disp,1);

O = [0,0];

tau = 0.5;

Bdispy = 0;

Cdispy = 0;

GG = zeros(2,N);

exitflags = zeros(1,N);

for i = 1:N

Peakdispx(i) = Disp(i,2);

Peakdispy(i) = Disp(i,3);

Adispy(i) = Peakdispy(i);

mslope(i) = (Disp(i,5)-O(2))/(Disp(i,4)-O(1));

Cdispx(i) = Peakdispx(i)-((Peakdispy(i)-Cdispy)/(mslope(i)));

Bdispx(i) = Cdispx(i)*-1;

Adispx(i) = (Adispy(i)-Bdispy)/(mslope(i))+Bdispx(i);

slopeACdisp(i) = ((Adispy(i)-Cdispy)/(Adispx(i)-Cdispx(i)));

Afull(i) = (0.5*abs(Adispx(i)*Bdispy+Bdispx(i)*Cdispy+Cdispx(i)*Peakdispy(i)+Peakdispx(i)*Adispy(i)-Adispy(i)*Bdispx(i)-Bdispy*Cdispx(i)-Cdispy*Peakdispx(i)-Peakdispy(i)*Adispx(i)));

gg0 = [tau*Disp(i,3);tau*Disp(i,2)];

[GG(:,i),exitflags(i)] = solveSys(Peakdispx(i),Peakdispy(i),Adispy(i),mslope(i),Cdispx(i),Cdispy,Bdispx(i),Bdispy,Adispx(i),slopeACdisp(i),Afull(i),tau,gg0);

%% pinch in x and y

disppinchy(i) = GG(2)/Disp(i,3);

U(i) = ((GG(2)-Cdispy)/mslope(i))+Cdispx(i);

L(i) = ((GG(2)-Bdispy)/mslope(i))+Bdispx(i);

disppinchx(i) = (GG(2)-L(i))/(U(i)-L(i));

end

function [GG,exitflag] = solveSys(Peakdispx,Peakdispy,Adispy,mslope,Cdispx,Cdispy,Bdispx,Bdispy,Adispx,slopeACdisp,Afull,tau,gg0)

[GG,~,exitflag] = fsolve(@fun,gg0);

function fval = fun(x)

% define function which will equal zero at desired values of x

% Note: your x = x(1) and your y = x(2)

fval(1) = -slopeACdisp+((gg0(2)-Cdispy)/(gg0(1)-Cdispx));

fval(2) = -tau+(0.5*abs((gg0(1)*Bdispy+Bdispx*Cdispy+Cdispx*Peakdispy+Peakdispx*gg0(2)-gg0(2)*Bdispx-Bdispy*Cdispx-Cdispy*Peakdispx-Peakdispy*gg0(1))))/Afull;

end

end

[Jon]

Sorry, I've haven't been on the computer lately. Did you find your problem yet?

Without digging into the details I see some things that could definitely cause problems for you.

Your function, defined at the bottom, fval = fun(x), does not acutally ever use the argument, x. The solver is going to try to vary the value of x to find where this results in fval(1) = 0 and fval(2) = 0. If there is no x in the expressions for fval(1) and fval(2) it can't possibly find a value of x that will cause these to be zero.

So you need to figure out what variable you would like the solver to vary in its search for fval=0 and include it.

All the other variables, slopeACdisp, Cdispy, etc are just parameters that are held constant while x varies. Also, your initial guess for the value of x, gg0 should not be included in the expressions for fval. fsolve uses this to assign the initial value for x in its hunt for the solution.

By the way, if you use the code button on the MATLAB answers toolbar you can cut and paste your code and have it nicely formatted.

Answer 2

[John D'Errico]

If these are LINEAR equations, then using fzolve or vpasolve is actively silly. They are designed to solve much more difficult problems, that is, nonlinear problems.

For something as simple as 2x2 matrices, you could even use inv. But backslash will be sufficient. For example:

N = 3;
a = [1,2,3];
b = [2,1,3];
d = [3,2,2];
e = [2,2,1];
C1 = [ 10, 4, 6];
C2 = [ 15, 3, 4];

M = reshape([a;d;b;e],[2,2,N]);
C = [C1;C2];

Now, you could trivially write a loop.

xy = zeros(2,N);
for n = 1:N
  xy(:,n) = M(:,:,n)\C(:,n);
end
xy

xy = 2×3

2.5000 2.5000 2.0000 3.7500 -1.0000 -0.0000

Could you do this without a loop? Well, yes. I could write it in a way that is fully vectorized. If N is not too large, then there is no need to introduce sparse matrices, but we could do so.

Mc = mat2cell(M,2,2,ones(1,N));
XY = blkdiag(Mc{:})\C(:);
XY = reshape(XY,2,N)

XY = 2×3

2.5000 2.5000 2.0000 3.7500 -1.0000 -0.0000

This latter solution did not use sparse matrices. But as I said, that is a waste of time for such a small problem.

Only if your problem is truly nonlinear would I bother with the other tools. And only if N is fairly large would I bother with sparse matrices, and even then only if you were solving this problem more than once. Note that backslash will factroizes the resultign block diagonal matrix into block tridiagonal matrices, so it will be fast. The only issue is if the memory required would become significant. Then the use of sparse tools would be a huge benefit. I'll not bother worrying unless you tell me your problem truly is linear, AND that N is significant.

[Jon]

Nice example of how you can solve this very cleanly if the equations are all linear.

In the OP's example the first equation was non-linear

a(i)*x^2+b(i)*y = C1(i);

due to the squared term. That is why I suggested using fsolve in my answer above