2 views (last 30 days)

Show older comments

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.

Walter Roberson
on 1 Jul 2021

syms Rab Rbc Rac x

Q = @(v) sym(v);

a = Q(10);

b = Q(30000);

gam111 = Gamcdf(x,a,b);

gam11 = 1-gam111;

a = Q(9);

b = Q(25000);

gam222 = Gamcdf(x,a,b);

gam22 = 1-gam222;

a = Q(12);

b = Q(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;

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

eq2

s2 = subs(eq2);

eq3

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

[S.Rab, S.Rbc, S.Rac]

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))

vpa(solve(gam2,x))

vpa(solve(gam3,x))

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)

Wait, is that a solution?

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

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)

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

Walter Roberson
on 5 Jul 2021

John D'Errico
on 1 Jul 2021

Edited: John D'Errico
on 1 Jul 2021

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)

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))

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.

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!