Solve system of equations and inequalities with multiple solutions?

How to solve a system of equations and inequalities with many solutions? I have the following set of equations:

a + b + c = 4 ;

9 <= a + 2b + 3c <= 12;

c >= 2;

b+c >= 3 ;

a >=0;

and get all the following solutions (a,b,c)=(1,1,2) or (0,2,2) or (1,0,3) or (0,1,3) or (0,0,4).

I have tried solving with vpasolve but it gives me "Empty sym: 0-by-1"

syms m1 m2 m3
S = vpasolve([m1 + m2 + m3 == 4, 9 <= m1+ 2* m2+ 3* m3, m1+ 2* m2+ 3* m3 <= 12, m3 >=2, m3+m2 >=3], [m1,m2,m3])

[Walter Roberson]

Also solutions in integers:

• -2 4 2
• -1 2 3
• -1 3 2
• 1 -1 4

and there are more solutions if non-integers are permitted, such as 0, 3/2, 5/2

[Hirak Basumatary]

@walter: Sir, i tried solving it with 'vpasolve' . I dont think i am trying the right technique. I need only integers. Also the value of a must be greater than equal to 0.

[Walter Roberson]

vpasolve() cannot deal with inequalities, usually. In cases that it does, it would only return a single solution, as inequalities are not polynomials.

You should be using

syms m1 m2 m3 integer
assumeAlso([m1 m2 m3] >= 0)

solve() can sometimes deal with inequalities. When it does work, it gives an "example" solution rather than an expression. But solve() cannot deal with that system in the way you have written it.

[madhan ravi]

linprog() ?

[Hirak Basumatary]

@madhan: linprog gives me the optimal solution. But how can i generate all the possible ositive real integer solutions

Answer 1

[Walter Roberson]

[M1, M2, M3] = ndgrid(0:4);
m1 = M1(:); m2 = M2(:); m3 = M3(:);
solidx = find(all([m1 + m2 + m3 == 4, 9 <= m1+ 2* m2+ 3* m3, m1+ 2* m2+ 3* m3 <= 12, m3 >=2, m3+m2 >=3], 2));
[m1(solidx), m2(solidx), m3(solidx)]

[Hirak Basumatary]

@walter: Sir, thank you very much. I will keep in mind your technique in the future.

Answer 2

[John D'Errico]

Remember that when you have inequalities, there are potentially infinitely many solutions. Even when you constrain the solutions to be integer, thus creating a Diophantine problem of some ilk, there may still be infinitely many solutions. Thus solving the problem

x - 2*y - z == 1

subject to positivity and integer constraints on x and y, you clearly cannot list all solutions. Admittedly, the solution to the simple one posed by me is pretty trivial. It looks like it could be expressed as

y = (1 - x + z)/2

which is valid only when x and z have opposite parity modulo 2, and x <= z + 1.

In general, solve is not designed to solve any fully general Diophantine problem. As for vpasolve, it is absolutely not designed to solve that class of problem. So if I were to pose some nasty class of Pell equation, I'd not expect it to provide the full set of solutions.

And no, tools like intlinprog are NOT a valid solution, even when the problem is fully linear, since as you observed, linear programming tools can only return one solution.

In the trivial example you pose, exhaustive enuneration is arguably best. Just generate all sets of integer triples that do not exceed 4, then test them all. Keep only those that satisfy all of your constraints. Walter shows you how that solution would apply here in his answer.

In a very large case however, this will be impossible. So if you have dozens of integer variables, each of which may be large, that exhaustive enuneration problem might involve far more possibilities than even a fast computer can handle.

Now, I do not know what your real target problem is here. People always pose trivially small problems as an example, when in their minds they are looking to solve something huge. Is your problem truly linear in all respects? If so, then there is much you want to learn in the field of linear Diophantine equations. Are you thinking of something like Pell equations, so now polynomial Diophantine equations of some general form? This will lead you down other back alleys in mathematics.

You can also look to use tools like recursion, if you are looking for a purely programmatic solution, but one that is still exhaustive.

For example, given the example problem, suppose we looked for all solutions where a == 0. Now your problem reduces to the simpler:

b + c = 4 ; 
9 <= 2b + 3c <= 12; 
c >= 2; 
b+c >= 3 ;

Do any solutions exist for b == 0?

c = 4 ; 
9 <= 3c <= 12; 
c >= 2; 
c >= 3 ;

We find that c == 4 is the only such solution.

Solutions = [0 0 4];

Return to the previous level, with that current solution. Can we increment b to b==1?

1 + c = 4 ; 
7 <= 3c <= 10; 
c >= 2;
c >= 2 ;

Now we find that only c==3 is a solution.

Solutions = [Solutions;0 1 3];

Increment b once more to b == 2?

c = 2 ; 
5 <= 3c <= 8; 
c >= 2; 
c >= 1 ;

Now c == 2 is the only solution.

Solutions = [Solutions;0 2 2];

Incrementing b once more to 3,

c = 1; 
3 <= 3c <= 6; 
c >= 2; 
c >= 0 ;

There we found an inconsistency, so no more solutions.

Next, increment a, then repeat the process. Each time, we reduce the problem, eliminating one variable. This process will resolve all possible solutions, as long as the set of solutions is finite, and not too large. Again, if you have hundreds of variables, with pentiillions of solutions, then you are up a creek with no paddle with any method you might hope to employ.

So, you may need to do some reading in the field of Diophantine eqautions, linear or nonlinear. Or you may be able to use recursive methods. For example, my partitions tool, found on the file exchange, is a very nice way to generate a list of all partitions of an integer. It was written to be recursive. But it will bog down for large problems, as you might expect. Or, your solution may simply be exhaustive enumeration. But the one thing that will usually fail if you try is to just throw it at solve or vpasolve. They simply are not designed to solve that class of problem.

Answer 3

[Basiya]

Solve and plot the system help me please

2x+6y+7z=8

6x+7y-3z<=5

-5x+9y03z>=7

x>=0

y>=0

z>=0

[Walter Roberson]

rng(123456)
syms x y z delta_2 delta_3

C = randi([-9 9], 3, 4)

C = 3×4

-7 8 -1 1 9 -2 6 -2 -5 -3 -7 -1

LHS = C(:,1:3)*[x;y;z]

LHS = 

RHS = C(:,4)

RHS = 3×1

1 -2 -1

eqns = [
    LHS(1) == RHS(1);
    LHS(2) + delta_2 == RHS(2);
    LHS(3) == RHS(3) + delta_3;
    ]

eqns = 

sol = solve(eqns, [x y z])

sol = struct with fields:
    x: (46*delta_3)/77 - (59*delta_2)/77 - 12/7
    y: (3*delta_3)/7 - (4*delta_2)/7 - 8/7
    z: (61*delta_2)/77 - (58*delta_3)/77 + 13/7

We can see from this that larger delta_3 implies larger x and larger y, but smaller z. Lower bound on delta_3 is

partial1 = solve(sol.x == 0, delta_3)

partial1 = 

partial2 = subs(sol, delta_3, partial1)

ans = struct with fields:
    x: 0
    y: 2/23 - delta_2/46
    z: - (4*delta_2)/23 - 7/23

and we can see from there that given that lower bound, that z will be negative unless delta_2 is negative. But to satisfy that >= inequality for the third equation, delta_2 and delta_3 must be non-negative. We conclude from this that the system of inequalities in this example has no solution.

Does the system of inequalities for your problem have a solution? Well, work it through.

Remember, expression1 <= expression2 implies 0 <= expression2 - expression1 which implies that there is some value non-negative delta such that delta == expression2 - expression1 which implies that expression1 + delta == expression2 . So you can transform the inequality into equality by putting in extra non-negative variables; and then you can later find the boundary conditions for those extra variables based upon the x>=0 and so on.