Asked by Tarmo Tukiainen
on 19 Nov 2019 at 9:41

I am trying to solve a cubic function from which I would like to get one numeric result as my variable. I've been able to do this by other means, but not with MATLAB and that's what I would need to do. I am not very good with Matlab for Im new with it so any help would be much appreciated!

My equation is; Solve for, V: V^3 - x * V^2 + y * V - z = 0.

In reality x = b+ ((R*T)/p), y = a/p and z = (a*b)/p. I know all these variables/results of the calculations so I only need to find out what V is.

So in one case the equation can be expressed as: V^3 - 0.10933127 * V^2 + 0.00082152 * V - 0.00002193 = 0

I have tried the following means and failed:

syms V

eqn = V^3 - 0.30993828*V^2 + 0.002598122*V - 0.0000693698 == 0;

S = solve(eqn,V,'Real',true)

With the first one I get symbolic results and cant open them by doubleclicking, I dont know if this is normal?

And

a = 1;

b = 0.30993828;

c = 0.002598122;

d = 0.0000693698;

function sols = solve_cubic(a, b, c, d)

syms x

sols = solve(a*x^3 - b*x^2 + c*x - d, x)

end

For some reason this does not give me anything for x in workspace, only values a,b,c,d that I already list.

And

p = [1 -0.10933127 0.00082152 -0.00002198];

r = roots(p)

With this I get a column vector with each row being "0.XXX + 0.XXXi".

If there is someone with extra knowledge thank you in advance!

Best regards,

-Tarmo

Answer by Walter Roberson
on 19 Nov 2019 at 10:02

Edited by Walter Roberson
on 20 Nov 2019 at 0:06

Accepted Answer

S = solve(eqn,V,'Real',true,'MaxDegree',3)

and

sols = solve(a*x^3 - b*x^2 + c*x - d, x, 'MaxDegree', 3)

and

p = [1 -0.10933127 0.00082152 -0.00002198];

r = roots(p);

r(imag(r) ~= 0) = [];

However, what you should be coding for the first equation is

syms V

syms d309 d259 d693 real

assume(-0.000000005 <= d309 & d309 < 0.000000005);

assume(-0.0000000005 <= d259 & d259 < 0.0000000005);

assume(-0.00000000005 <= d693 & d693 <= 0.00000000005);

eqn = V^3 - (0.30993828+d309)*V^2 + (0.002598122+d259)*V - (0.0000693698+d693) == 0;

S = solve(eqn,V,'Real', true,'MaxDegree', 3, 'returnconditions', true);

Now S.V will be the generalized solution for V, and S.conditions will be the conditions under which the solution is valid. Both of them will be in terms of the variables expressing the uncertainty in your floating point values, since (for example) 0.30993828 does not mean 30993828/100000000 exactly in an equation: in an equation, 0.30993828 expresses "the set of valeus that rounds to 30993828/100000000 exactly". Depending on the exact difference between the set member and the nominal representative 30993828/100000000, the real root might not exist !

If you have a floating point coefficient, then you probably should not be using solve(). solve() is reserved for finding indefinitely precise closed-form solutions whenever possible. If your desired solution is double precision numbers such as 0.302098127847416 then you should not be using solve(). Use roots() or use vpasolve() instead.

Tarmo Tukiainen
on 19 Nov 2019 at 12:22

Thank you both for your anwers, what solved it was that one line of code you gave me, Walter.

r(imag(r) ~= 0) = [];

Now for my r value I get the same one positive numeric value as given to me by different a calculator/equation solver for V.

I tried it with different values for temperature and pressure and I get the same values as on paper.

So now my code looks like this in case someone needs to solve a similar problem later.

% Input data(a,b,R,Temp,p) located above in script

A = 1 ;

B = -(b + ((R*Temp)/p)) ;

C = a/p ;

D = -((a*b)/p) ;

%% Solving equation

% A*V^3 - B*V^2 + C*V - D = 0

p = [A B C D]

r = roots(p);

r(imag(r) ~= 0) = [];

Sign in to comment.

Answer by ME
on 19 Nov 2019 at 10:01

Edited by ME
on 19 Nov 2019 at 10:03

You aren't going to solve this properly by getting an output of one numerical solution because a cubic equation does actually have three solutions.

Of your options I would use the roots function. My hunch is that if those aren't the correct solutions to the cubic then you have probably not defined your coefficients correctly before putting them into roots.

If you are after a single solution then you may need to take the output and apply some other functions to pick the one you want. For example if you wanted the real valued solution only then you could use:

r=real(roots(p))

or something to that effect.

Sign in to comment.

Opportunities for recent engineering grads.

Apply Today
## 0 Comments

Sign in to comment.