how to solve 3 equations with a known variable and 3 unknown variables

I want to solve this equation in Matlab:

m=0:2000;

(2-x)+(3-y)z=f1(m)

xy−3z=f2(m)

x−(2-y)(5-z)=f3(m)

I tried with solve but I didnt get the correct result. Could you please help me with the code?

And then I want to plot them.

[John D'Errico]

Why not show what you tried? If you don't, then we need to do all of your work for you. Start by defining m as a parameter, thus

syms x y z m

Do NOT define m in advance as a vector. After you have solved the system for x,y,z, paramerically as a function of m, only then should you substitute the sequence of values for m.

So try it.

[farid p]

Could you please read the edited question again?

[John D'Errico]

Should I keep on chasing a moving target? It seems you have now created some fully generic problem, with functions of m in there where you will not even admit what they are.

You have a nonlinear system of equations, in three variables, where the right hand side is a function of m. So?

First, write down the functions f1(m), f2(m), f3(m). If you cannot express tham then in some functional form, then what can you do? Nonlinear systems of equations are not assured to always even have a solution.

SOLVE the system, for x,y,z as a function of the parameter m. Note that since your equations are nonlinear, you may easily have bifurcations, all sorts of nasty things happening. But assuming a solution does exist...

Finally, substitute m into the solution. Now you can plot the result, AS A FUNCTION OF m.

So what is the problem? The solution is exactly as I described it before, though since the system is nonlinear, all sorts of nasty stuff may happen. As I said, bifurcations, etc.

Will you now change the question again? TRY IT. Show what you tried.

[farid p]

%Hi.
%here is my exact code; however it gives an error, because gamcfd works only with variable of x and I can not provide %it with " sym x". Could you please help me again.

clc
clear all
close all

syms Rab Rbc Rac x
a = 10;
b = 30000;
gam111 = gamcdf(x,a,b);
gam11 = 1-gam111;
a = 9;
b = 25000;
gam222 = gamcdf(x,a,b);
gam22 =  1-gam222;
a = 12;
b = 45000;
gam333 = gamcdf(x,a,b);
gam33 = 1-gam333;
gam1=gam11.*gam22;
gam2=gam22.*gam33;
gam3=gam11.*gam33;
eq11=1- ( (1-Rab)*(1-(Rbc*Rac)))==1; 
eq22=1- ( (1-Rbc)*(1-(Rab*Rac)))==1; 
eq33=1- ( (1-Rac)*(1-(Rbc*Rab)))==1; 
eq1=eq11.*gam1;
eq2=eq22.*gam2;
eq3=eq33.*gam3;
x=0:56000;
eq1
subs(eq1)
x=0:56000;
eq2
subs(eq2)
x=0:56000;
eq3
subs(eq3)
S=solve(eq1,eq2,eq3,Rab,Rbc,Rac);
plot(x,s)

[farid p]

clc
clear all
close all

syms Rab Rbc Rac x
a = 10;
b = 30000;
gam111 = gamcdf(x,a,b);
gam11 = 1-gam111;
a = 9;
b = 25000;
gam222 = gamcdf(x,a,b);
gam22 =  1-gam222;
a = 12;
b = 45000;
gam333 = gamcdf(x,a,b);
gam33 = 1-gam333;
gam1=gam11.*gam22;
gam2=gam22.*gam33;
gam3=gam11.*gam33;
eq11=1- ( (1-Rab)*(1-(Rbc*Rac)))==1; 
eq22=1- ( (1-Rbc)*(1-(Rab*Rac)))==1; 
eq33=1- ( (1-Rac)*(1-(Rbc*Rab)))==1; 
eq1=eq11.*gam1;
eq2=eq22.*gam2;
eq3=eq33.*gam3;
x=0:56000;
eq1
subs(eq1)
x=0:56000;
eq2
subs(eq2)
x=0:56000;
eq3
subs(eq3)
S=solve(eq1,eq2,eq3,Rab,Rbc,Rac);
plot(x,s)

Answer 1

[Walter Roberson]

syms Rab Rbc Rac x
Q = @(v) sym(v);
a = Q(10);
b = Q(30000);
gam111 = Gamcdf(x,a,b);

p = 

gam11 = 1-gam111;
a = Q(9);
b = Q(25000);
gam222 = Gamcdf(x,a,b);

p = 

gam22 =  1-gam222;
a = Q(12);
b = Q(45000);
gam333 = Gamcdf(x,a,b);

p = 

gam33 = 1-gam333;
gam1=gam11.*gam22;
gam2=gam22.*gam33;
gam3=gam11.*gam33;
eq11=1- ( (1-Rab)*(1-(Rbc*Rac)))==1; 
eq22=1- ( (1-Rbc)*(1-(Rab*Rac)))==1; 
eq33=1- ( (1-Rac)*(1-(Rbc*Rab)))==1; 
eq1=eq11.*gam1;
eq2=eq22.*gam2;
eq3=eq33.*gam3;

Let us look at those equations for a moment. eq11, eq22, eq33 are independent of x.

All three of them are equations, and you multiply the equations by expressions that involve x. The multiplication acts on both sides of the equation, so if you had f(x) == c and multiply both sides by constant A then A*f(x) == A*c would be the resulting equation.

The multiplication does not introduce new solutions: the resulting extended equation has solutions either when the original equations (eq11, eq22, eq33) are true by themselves, or when the multipliers (gam1, gam2, gam3) are 0. So because eq11, eq22, eq33 are independent of x, the multiplied values, eq1, eq2, eq3, will have solutions at the points where solve([eq11, eq22, eq33], [Rab, Rbc, Rac]) works, and will also have solutions where solve([gam1, gam2, gam3], x) works.

%x=0:56000;
eq1

eq1 = 

eq2

eq2 = 

s2 = subs(eq2);
eq3

eq3 = 

S = solve(eq1,eq2,eq3,Rab,Rbc,Rac);
[S.Rab, S.Rbc, S.Rac]

ans = 

And there you see the solutions that are independent of x, as discussed above. You could plot them over x, but it would be constant graphs.

Okay, how about the case where the x contributions become 0, solve([gam1, gam2, gam3],x) ?

vpa(solve(gam1,x))

ans = 

vpa(solve(gam2,x))

ans = 

vpa(solve(gam3,x))

ans = 

You can see that for each individual gam* that most of the solutions are imaginary, and all of them are outside x in 0:56000 . You can also see that although the first of them shares a real solution with the third, and the second of them shares a real solution with the third, that there is no real-valued x for which all three have a solution. Therefore when you solve the three of them as a group, you are not going to see any real-valued solutions.

... Not unless solve() has omitted some solution. So let us try:

G = sum([gam1,gam2,gam3].^2);
Gsol = vpasolve(G);
vpa(Gsol,16)

ans = 

1.624822280811449e+25

Wait, is that a solution?

subs([gam1,gam2,gam3],x,Gsol)

ans = 

Looks like it??

... But no, it is a false solution, where the values reach 0 because of round-off, but are non-zero. See

digits(1000)
Gsol2 = vpasolve(simplify(G));
vpa(Gsol2,16)

ans = 

8.875782315672845e+1014

More digits, the solution extends out a lot further... in fact there is no solution until you reach x = infinity.

function p = Gamcdf(x,a,b)
   p = Gammainc(x/b, a)
end

function y = Gammainc(x,a)
   syms t real
   y = int(t^(a-1)*exp(-t), t, 0, x) / gamma(a);
end

[farid p]

Yeah i actually want Rab, Rbc, Rac if they are between ero and one; like:

0< Rab < 1

0< Rbc < 1

0< Rac < 1

Also, they should be real numbers (not complex number).

[farid p]

Oh yes. this is the code that I wrote in the first place:

[farid p]

a = 7;
b = 7*(10^4);
syms Rab Rca Rbc x
y_gam11 = gamcdf(x,a,b);
y_gam1=1-y_gam11;
a = 0.1*7;
b = 1.1*5.45*(10^4);
y_gam22 = gamcdf(x,a,b);
y_gam2=1-y_gam22;
a = 0.075*7.06;
b = 1.3*5*(10^4);

y_gam33 = gamcdf(x,a,b);
y_gam3=1-y_gam33;

gammaab=y_gam1.*y_gam2;
gammabc=y_gam2.*y_gam3;
gammaca=y_gam3.*y_gam1;
eq1=( 1-( (1-Rab)*(1-(Rca*Rbc)) ) )==gammaab;
eq2=( 1-( (1-Rbc)*(1-(Rab*Rca)) ) )==gammabc;
eq3=( 1-( (1-Rca)*(1-(Rab*Rbc)) ) )==gammaca;
s=solve(eq1,eq2,eq3,Rab,Rbc,Rca);

vpa(s,Rab)
vpa(s,Rbc)
vpa(s,Rca)
clc
clear all
close all
a = 7;
b = 7*(10^4);
syms Rab Rca Rbc x
y_gam11 = gamcdf(x,a,b);
y_gam1=1-y_gam11;
a = 0.1*7;
b = 1.1*5.45*(10^4);
y_gam22 = gamcdf(x,a,b);
y_gam2=1-y_gam22;
a = 0.075*7.06;
b = 1.3*5*(10^4);

y_gam33 = gamcdf(x,a,b);
y_gam3=1-y_gam33;

gammaab=y_gam1.*y_gam2;
gammabc=y_gam2.*y_gam3;
gammaca=y_gam3.*y_gam1;
eq1=( 1-( (1-Rab)*(1-(Rca*Rbc)) ) )==gammaab;
eq2=( 1-( (1-Rbc)*(1-(Rab*Rca)) ) )==gammabc;
eq3=( 1-( (1-Rca)*(1-(Rab*Rbc)) ) )==gammaca;
s=solve(eq1,eq2,eq3,Rab,Rbc,Rca);

vpa(s,Rab)
vpa(s,Rbc)
vpa(s,Rca)

[Walter Roberson]

Your latest version of the equations has no solution in (0,1) .

You can take

eq = [eq1, eq2, eq3]
err = sum((lhs(eq)-rhs(eq)).^2);
err10 = subs(err,x,1/10)
ERR1 = matlabFunction(err10, 'vars', {[Rab, Rbc, Rca]})

and if you now ask to fsolve(ERR1) you might get results. However if you look at the output at the locations given, it is non-zero, and that is not just a matter of round-off.

If you proceed

s1 = solve(diff(err10,Rab),Rab);
e2 = subs(err10,Rab,s1);
s2 = solve(diff(e2,Rbc),Rbc);
e3 = subs(e2,Rbc,s2);

At this point, e3 will be a vector of 5 results involving root() of polynomials of degree 5.

If you then try to

s3 = solve(diff(e3(1),Rca),Rca)

it will not succeed.

When I look at the e3 values in Maple, I can see that e3 decreases as Rca goes from 0 to 1, until it reaches about 0.999058 at which point there is a discontinuity that rises abruptly between 0.999058 and 0.999059 . The e3 values are small, down to 2E-22 near 0.9990583, but they do not seem to vanish (though at the moment I cannot rule out the possibility that they get indefinitely close to 0 if you go indefinitely precise)

[farid p]

How can I define variables, so the results converge?

Thank you for your support.

[Walter Roberson]

I do not understand your question? Your equations are whatever they are, and what they are is not solvable. In order to have something solvable you would need different equations, or at least different constants. From where I am, I have no idea which of your constants are fixed and which of them can be adjusted.

Answer 2

[John D'Errico]

The solution is simple, as I said, and somehow, I imagine you will end up changing the question again. I'll give you an answer, to the question as it is now (version 2). However, I will not bother to solve the 3 variable problem, as it is just a more complicated version of a 1 variable problem. And since you seem to be just changing the equations arbitrarily, this means your question is just exemplar.

For example, suppose you want to solve the problem

x^2 + x = sin(m)

First, formulate the problem in a symbolic form. I DID TELL YOU TO DO EXACTLY THIS.

syms x m
eqn = x.^2 + x == sin(m);
xsol = solve(eqn,x,'maxdegree',2)

xsol = 

Because the problem is polynomial in x, there are two solutions, but you need to recognize that those roots may turn complex. If you care only about real solutions, they essentially disappear. But if your problem is fully nonlinear, it may be impossible to find an analytical solution. Since you seem insistant that the equations may be fully general, then you will need to learn about solutions of nonlinear equations, and how to solve them. Using a numerical solver can easily have problematic results. You will need to understand concepts like the basins of attraction for a nonlinear problem.

But, for the above problem, assume that m varies on the interval [0,2*pi].

fplot(xsol(1),[0,2*pi])

I could have added a plot for xsol(2) on the same axes. Note that the solution dissppears on the interval roughly [3.3,6] for m. That is because the solution becomes complex on that sub-interval. Thus,

vpa(subs(xsol(1),m,4))

ans = 

For much nastier general equations, expect that sort of thing to happen, but remember that a nonlinear equation may not have an algebraic solution in any nice form. Consider that even a general 5th degree polynomial has no such solution, and that a system of polynomial equations will often be equivalent to a higher degree polynomial equation.

Since you want to know how to solve this in the fully general case. I cannot give you a magical answer that will always work, since the problem may be arbitrarily nasty.

[farid p]

clc
clear all
close all

syms Rab Rbc Rac x
a = 10;
b = 30000;
gam111 = gamcdf(x,a,b);
gam11 = 1-gam111;
a = 9;
b = 25000;
gam222 = gamcdf(x,a,b);
gam22 =  1-gam222;
a = 12;
b = 45000;
gam333 = gamcdf(x,a,b);
gam33 = 1-gam333;
gam1=gam11.*gam22;
gam2=gam22.*gam33;
gam3=gam11.*gam33;
eq11=1- ( (1-Rab)*(1-(Rbc*Rac)))==1; 
eq22=1- ( (1-Rbc)*(1-(Rab*Rac)))==1; 
eq33=1- ( (1-Rac)*(1-(Rbc*Rab)))==1; 
eq1=eq11.*gam1;
eq2=eq22.*gam2;
eq3=eq33.*gam3;
x=0:56000;
eq1
subs(eq1)
x=0:56000;
eq2
subs(eq2)
x=0:56000;
eq3
subs(eq3)
S=solve(eq1,eq2,eq3,Rab,Rbc,Rac);
plot(x,s)