Which routine to be used for numerically solving a three-order equation represented in an array?

Hi all,

I'm sorry if the title is not clear enough. What I want to solve is an equation like

y^3 + a y^2 + b y + c = 0,

where y=y(x) and c is an array (or function) of x (I don't have an analytical expression for it) and for a,b they are constants. I know the fsolve routine, but as far as I understand, it can only be used when the explicit expressions of the coefficients are known. How to solve this polynomial? Thanks a lot.

[Walter Roberson]

I am not clear what is to be solved for? Is it the case that you know y and you know c, and you want to find a and b?

[Mengqi]

Hi, I don't know y, I want to solve y. c is an array actually, so I want to have a numerical result of y. a and b are scalar constants.

Answer 1

[Torsten]

For x fixed, you know a, b and c.

So for this fixed x, determine y(x) as

y=roots([1 a b c])

Best wishes

Torsten.

[Mengqi]

Sorry, roots seems not work. c is an array actually, and I want to have a numerical result of y. a and b are scalar constants. If c is a scalar, roots works.

Answer 2

[John D'Errico]

Trivially, one could just use a loop! Loops are not the worst thing in the world.

No matter what you do, the problem is that you will get three solutions for every element of c. So you will need to choose the root of interest.

Or you could use arrayfun.

For example...

a = 2;
b = 5;
c = magic(3);
res = arrayfun(@(ci) roots([1 a b ci]),c,'uniformOutput',false);
res
res = 
  [3x1 double]    [3x1 double]    [3x1 double]
  [3x1 double]    [3x1 double]    [3x1 double]
  [3x1 double]    [3x1 double]    [3x1 double]

res{1,1}
ans =
       -0.123914110557906 +      2.13316845986304i
       -0.123914110557906 -      2.13316845986304i
        -1.75217177888419 +                     0i