Solve a cubic equation using MATLAB code

I have a cubic equation whose coefficients are varying according to a parameter say w in the following manner:

a=2/w;
b=(3/w+3);
c=(4/(w-9))^3;
d=(5/(w+6))^2;

a*(x^3)+b*(x^2)+c*x+d=0

I want to solve the above equation using a m-file not in the command window. Is it possible ? I am fairly new to MATLAB, so help would be appreciated.

Answer 1

[Matt Fig]

From your comments, it looks like you want this instead:

function [] = solve_cubic()

w=10:2:30;
a=4*((9*w)-7);
b=(-2)*((18*w)-16).*(w+1);
c=3*((w+1).^2).*((3*w)-1);
d=(-3)*((1+w).^3);
f = zeros(1,length(w));

for ii = 1:numel(w)
    g = roots([a(ii),b(ii),c(ii),d(ii)]);
    f(ii) = g(g<1);
end

ptr=2./(1+w);
u=(1+w)./2;
l=(u.*(1-ptr.*f).^3)./3;
plot(w,l);

[Bhagat]

Why did you create the function solve_cubic ? It serves no purpose here,does it ?

[Bhagat]

w=10:2:30;
a=4*((9*w)-7);
b=(-2)*((18*w)-16).*(w+1);
c=3*((w+1).^2).*((3*w)-1);
d=(-3)*((1+w).^3);
f = roots([a,b,c,d]);
f=f(f<1);
ptr=2./(1+w);
u=(1+w)./2;
l=(u.*(1-(ptr.*f).^3))./3;
plot(w,l);

============================

In the above code following error is coming:

??? Error using ==> times

Matrix dimensions must agree.

Error in ==> cubic at 10

l=(u.*(1-(ptr.*f).^3))./3;

==================================

How and which matrix dimensions don't agree here ???

[Walter Roberson]

Your a, b, c, and d are all row vectors. [a,b,c,d] is then a row vector of size equal to the sum of the sizes of a, b, c, d. Your w vector has 11 elements, so your a, b, c, and d are all 11 long, so the constructed vector is 11+11+11+11 = 44 elements long. And then roots() will calmly tell you the 43 roots of that 44'th degree polynomial.

I told you above that roots() must be applied to each [a(K),b(K),c(K),d(K)] combination individually.

Read the documentation for roots(): it clearly says that the input must be a row vector. No possibility of a matrix of rows to be solved independently is mentioned.

If you choose to believe that roots() might accept a matrix, then go ahead and try

f = roots([a(:),b(:),c(:),d(:)])

as that will at least have one row per value of w.

Also: You are throwing away some of the solutions for f. How many solutions will you have left? If for each w, you do not happen to have exactly one real-valued root with value less than 1, then the size of f will not match the size of the ptr matrix that you are trying to multiply by.

In order to use your formula for "l", for each value of f, you must use the value of w that was used to generate that f. I warned you above about having 0, 1, or 3 roots per w...

[Matt Fig]

You keep posting the same erroneous code, but did you try to copy what I wrote and see if it gives you the values you expect? I ran it and it got no errors. If you don't want it as a function, simply delete the first line; it isn't hard to do!

[Ubong]

Seen the final response by Matt, that was just so nice . Still have questions for Matt though

1. what is the consequence of not adding the f = zeros(1,length(w)); line? I found that I can still get my row vector, f, without adding it

2. What do I do, if g had returned both real and complex numbers for every ii, and I need to get only real numbers from g

Answer 2

[Bruno Luong]

Please download this FEX,

There is a file called CardanRoots.m

w = linspace(-3,3,10);
a = 2./w;
b = (3./w+3);
c = (4./(w-9)).^3;
d = (5./(w+6)).^2;

X = CardanRoots([a(:) b(:) c(:) d(:)])

Each row of X contains three solutions of the respective cubic equation. Filter out the complex solutions if you don't need them.

[Ann Ytterberg]

Where is this CardonRoots.m?

[Image Analyst]

It's the first hit from Google: CardanRoots

[Bruno Luong]

Click on the link it should work

Answer 3

[Walter Roberson]

EDIT: corrected missing syms, added more detail

function sols = solve_cubic(a, b, c, d)
  syms x
  sols = solve(a*x^3 + b*x^2 + c*x + d);
end

The inputs to this can include symbolic expressions.

The outputs of this will be symbolic numeric radicals if the inputs are all numeric, but might include symbolic RootOf expressions if any of the arguments are symbolic. To convert symbolic numbers to double precision floating point numbers, use double() on the output.

[Bhagat]

I wrote the following code:

function f = e(a, b, c, d);
w=10:2:30;
a=2./(w+1);b=3./(w+2);c=4./(w+4);d=20;
f = solve(a*(x^3) + b*(x^2) + c*x + d);

---------------------------------------------

It says error 'need to define x'. Why would i define x when it is the variable i need to find ? And if i have to then how ?

[Bhagat]

ok ok i got it. the expression was not in inverted commas.

But the output this is giving in still in terms of a, b ,c,d. I want the actual value. Which function to use ???

[Walter Roberson]

Sorry I missed a statement before.

[Walter Roberson]

Don't use the apostrophes, use the syms command.

[Bhagat]

ok, but i got the values using apostrophes also. reason for using syms command

Answer 4

[Walter Roberson]

Provided that the input values are scalar numeric:

function sols = solve_cubic(a, b, c, d)
  sols = roots([a,b,c,d]);
end

EDIT: provide more detail

The solutions to this will be a column vector of double precision floating point numbers (possibly complex.)

[Bhagat]

Thanks i got the values. I need a bit more help though.

Can I use 'if' statement here so as to plot the values of the output greater than 1 with the parameter value 'W' ? how ?

[Walter Roberson]

Wvals = 1:exp(1):pi^3; %chose appropriate vals
nW = length(Wvals);
out1 = nan(1,nW);
out2 = nan(1,nW);
out3 = nan(1,nW);
for K = 1 : nW
 w = Wvals(K);
 a=2/w; b=(3/w+3); c=(4/(w-9))^3; d=(5/(w+6))^2;
 these_sols = solve_cubic(a,b,c,d);
 these_sols = these_sols(imag(these_sols) == 0);
 L = length(these_sols);
 if L > 0
 out1(K) = these_sols(1);
 end
 if L > 1
 out2(K) = these_sols(2);
 end
 if L > 2
 out3(K) = these_sols(3);
 end
end
plot(Wvals, out1, '.', Wvals, out2, '.', Wvals, out3, '.');

This will plot up to 3 points per w value, as there can be up to 3 real solutions in general. (I would have to look further to see whether 3 real solutions are actually possible.)

[Walter Roberson]

There is an algebraic theorem that any cubic in real coefficients has either one or three real roots, never 0 or 2. You can use that theorem to simplify the above code slightly.

I tried some values myself, and found that indeed for most values of w, there is only one real root. There is a region, though, that is somewhere within -10 to +10, in which there are 3 real roots for each value of w.

[Bhagat]

Ok, thank you. Now,Is there any function to find the values in a matrix which is less than a particular no. (not using the if statement)???

[Walter Roberson]

find(Matrix < ParticularNumber)

[Bhagat]

find(matrix<particular number) not working. It is giving the values 4 to 43 in steps of 1 represented in matrix form.

[Bhagat]

find statement gives the index of the value following the condition,not the value itself. Is there a statement to find the value less than a particular no. or ??

[Walter Roberson]

Matrix(Matrix < ParticularNumber)

[Bhagat]

Thx a lot. Will run my program and let you know if have any other doubts. You have been great

[Bhagat]

I wrote the following code:

------------------------

w=10:2:30;
a=4*((9*w)-7);
b=(-2)*((18*w)-16)*(w+1);
c=3*((w+1)^2)*((3*w)-1);
d=(-3)*((1+w)^3);
f = roots([a,b,c,d]);
f=f(f<1);
ptr=2./(1+w);
u=(1+w)./2;
l=(u*(1-ptr*f)^3)./3;
plot(w,l);

--------------------------

It shows the following error:

??? Error using ==> mtimes

Inner matrix dimensions must agree.

Error in ==> cubic at 3

b=(-2)*((18*w)-16)*(w+1);

======================

I am not able to correct it .Please help........

[Bhagat]

If i put w as a constant value, then it works but not it shows the error mentioned when I vary w.

[Walter Roberson]

You _are_ putting in w as a constant value: it's value is the entire list 10:2:30 . Using a list like that is not varying w.

Although it is possible to slightly change your a, b, c, and d so that values for each of the w would be calculated, the roots() call will not accept a vector of values, and must be called for each individual a, b, c, d combination.

w0 = 10:2:30;
for K = 1 : length(w0)
 w = w0(K);
 %calculate a, b, c, d here, and calculate the roots and select those in the proper range and continue on through to the calculation of "l". Do _not_ do the plot() at this time though.
 L{K} = l;
end

Now you can plot of w0(K) against L{K}, but keep in mind that you might have no real roots or one real root or three real roots for any individual w, so you will have to build appropriate plotting code.

[Bhagat]

Can't we accept more than 1 answer ?

Answer 5

[Awa Kologo]

1. Solving for the roots of a third order polynomial (i.e. finding a value that will make the expression equal zero) may require tedious algebra to do by hand but can be solved easily by a computer using an iterative approach. For example, finding the roots of the expression: , (ie. a value of x so that the equation is satisfied) is time consuming to do by hand. However, plugging in a guess for and then modifying that guess until a tolerance is met gives . Write a MATLAB script that solves exponential equations of the form where are constants that the user inputs in a 1 x 4 row vector: [a b c d]. The program will first check whether or not the input is a 1 x 4 matrix of numerical entries throw an error if it is not. The program will solve for one root of the polynomial iteratively to 6 decimal places and print out the value. You can take this to mean that your program should keep running for as long as Your program will need to make an initial guess for in order to solve iteratively. In most cases, the will give a reasonable first guess. You can use the roots() function to check your answer.

[Walter Roberson]

This does not appear to be an Answer to the Question that was asked.