Depressed cubic polynomial

Function to obtain the "depressed" cubic polynomial, from the general cubic form.

Description

Syntax:
   p2 = depressedcubic(p1)

If the general cubic polynomial p1(y) has the form p1=[a b c d], meaning:

 p1(y) = a*y^3 + b*y^2 + c*y +d

The Matlab function depressedcubic.m transform it into the "depressed" cubic form p2(x)

p2(x) = a*x^3 + k1*x + k2

by applying the change of variable

y = x - b/(3a)

where k1 = c - b^2/(3a)
and k2 = d - bc/(3a) + (2*b^3)/(27a^2)

This function checks if the input is a valid cubic polynomial (expressed as a 4-vector as in Matlab format) with the coefficient of the cubic term non zero.

The method is the work of Nicolo' Fontana (a.k.a. Tartaglia; 1500-1557) and Scipione Del Ferro (1465-1526). Later published by Cardano in his Ars Magna.

• Required Products: MATLAB
• MATLAB release: MATLAB 7.1.0 (R14SP3)

[Isaac Mancero Mosquera]

"Thank you for the feedback. I will improve this script with a better help section and possibly a change in the name" but I'm repeating myself here.

I don't feel attacked if you criticize my program, I'm actually thankful, since programming is not my cup of tea. Also, I always welcome critics; working in Science, I have to be prepared to debunk and be debunked, there is nothing personal in it; but if we disagree in something, we have to proceed separate ways -still nothing personal in it.

I'm uploading a new version of the file soon.

[John D'Errico]

I agree. Trying to help you when you refuse to learn, trying to show you how to make your code better is a silly thing when you refuse to do so. I should have known better.

[Isaac Mancero Mosquera]

Attacking a mathematically silly idea is not attacking you. Stop the drama.

[John D'Errico]

Sigh. Clearly you have not read a word of what I said, while at the same time stating that I did not read your response. Is this one of those attempts to bamboozle the reader by attacking me, claiming I have done what you are doing yourself?

I said that it makes sense to offer TWO options. One that works as you like and one that works as others will surely want. If you bothered looking at the pseudo-code I offered, you will see it does exactly that. This way the user can have it either way. You offer only ONE option, which makes your code less useful than it might be. Flexible software is useful software.

Sorry, but I see no reason to waste any more of my time on this, as you already know the truth.

[Isaac Mancero Mosquera]

I'll try to be more clear, for the sake of Math (I agree with all your observations regarding the programming). This two polynomials:

p(x) and q(x)=K*p(x) are different.

They even have the same roots:

p(x)=0 and q(x)=0 at the same "x".

But they are not the same as functionals, and this is best seen when wee equal them to "w" non zero:

p(x)=w & q(x)=w have not the same solutions.

Multiplying polynomials by constants is a very popular method when finding roots, or when you have the other side of the equation p(x)=w, then:

K*p(x) = K*w , still mathematically correct.

In my program, there is no a hint of any equation, is just a mapping of a polynomial from one form to other.

That's why I put a warning to the user, they might want to multiply and normalize, but they do that themselves, at their own risk; because this program does not have "w" (from p(x)=w) as an input.

[John D'Errico]

I suppose I have no choice but to disagree with you.

[Isaac Mancero Mosquera]

Sure. If that coefficient is not 1, you can't apply this method and the code is not useful. What's so wrong with that? Are you aware that Not all theorems apply when you want, but when their logic premises hold true?

If the purpose of this file was the mapping of general polynomial forms, then you'd be right and this would be a different story.

You are not reading my answer well: The depressed form can be used in a NON-EQUATION context, meaning not in a problem of the form: p(x)=0. In THAT case, it's better for the user to normalize whatever thing they want by themselves.

Please note:

Roots of p(x)=0 are unchanged by multiplication by a scalar.

Roots are changed in the case of p(x)=K. If you normalize p(x), then you should normalize K. Since K is nowhere to be found in my file, that task should be left to the user.
[Because delferro.m has a more general mapping of a polynomial unrelated to root-findings methods, so you shouldn't go around happily multiplying polynomials!].

You may like to multiply coefficients just because you want. But you are not supposed to do so in whatever context.

IF (big IF) your problem is related to findings roots, then it would be ok, but then again, there is already a built-in Matlab function to do that.

IF (again big IF) you are dealing with the functional w=p(x), you can normalize BOTH sides in your equation, or apply a more general solution, which is not exactly "divide by the coefficient of x^3", but applying another change of variable to obtain the monic trinomial form.

I hope it's clear now why this function will not divide by the coefficient of x^3 (the general mapping formulation won't do it either).

[John D'Errico]

Essentially you are saying that since you sometimes use this tool in a case where the first coefficient is NOT unity, you prefer the way it currently behaves. That is silly in the extreme, since the code returns NOTHING (but a flag that is not even a warning, and is certainly not an error) if the first coefficient is not unity. This code is useless to even you if that first coefficient is not unity. So put in the documentation that it REQUIRES a unit first term. Or put a flag in that allows the user to change the behavior of the code as they need it.

For example, do this...

================================
function Pd = depressed(P,flag)
% H1 line and help, that fully documents the arguments and any implicit assumptions.

% do we need to provide a default for flag?
% note that defaults must be stated in the help.
if (nargin < 2) || isempty(flag)
flag = true;
end

% assure that P(1) is unity.
if flag
% scale the poly, to be later undone
p1 = P(1);
P = P/p1;
else
% the user has chosen to enforce a unit first
% coefficient. If they forgot to do that, but still
% called this code, get upset and refuse to work.
if P(1) ~= 1
error('The sky is falling')
else
p1 = 1;
end
end

% Compute the Del Ferro depressed form
k1 = (stuff);
k2 = (more stuff);

% Build the depressed polynomial, and undo the scaling of P.
% If flag was false and we escape the error condition then
% p1 == 1 and no scaling was done, so the multiply is
% an irrelevant sham.
Pd = p1*[1 0 k1 k2];
================================

Note the internal comments that explain what I'm doing in each block of code. Those comments are crucial for yourself and for others. They help when you need to debug the code in the future. They help the user to understand what you are doing. Get used to writing these comments even as you write code, as great memory aids to yourself. Make it an unconscious habit.

[Isaac Mancero Mosquera]

Thank you for the feedback. I will improve this script with a better help section and possibly a change in the name.
-Dividing for the first coefficient after testing that it's not zero- is fine if you're going to find roots in p(x)=0. The depressed cubic for is interesting in other contexts too.
In my case, the need to invert a cubic polynomial regression: y=p(x), such that the roots of p(x)-y=0 are needed, and I prefer the user to be aware of the need of a normalization (since it will also divide potential user data "y")
There is however another more general formula, I'll include in a future version.

[John D'Errico]

 

[John D'Errico]

There are a few problems here, not terribly bad problems, but worth correcting.

First of all there is no H1 line. While there is a help block, the first line of the help should be a descriptive line, that includes some useful keywords. It should help the user find this tool when they have no idea next month what the name of the function is. Perhaps this:

% delferro - convert cubic polynomial to a depressed form

See that this includes a few well chosen keywords, so that when you cannot remember delferro as the name, you have a remote chance of finding this tool with the function lookfor. (Learn to use lookfor.)

Ok, next, I'm not at all in love with the name delferro. The name of a function is best given as a descriptive word, that will help you remember what it does. I might suggest depressedform, or just depressed.

You might also provide a link, perhaps this one for documentation: http://en.wikipedia.org/wiki/Cubic_function

How about the code? It is pretty short, so here it is:

-------------------------------------
if ( p1(1) ==1 )
a = p1(2);
b = p1(3);
c = p1(4);

k1 = b - (1/3)*(a.*a);
k2 = c - (1/3)*(a.*b) + (2/27)*(a.*a);

p2 = [1 0 k1 k2];

else
disp('ERROR: coefficient "a3" of x^3 should be 1')
disp('Advice: you can enter p(x) / a3')
p2 = [];
end;
-------------------------------------

(No need by the way to have a semi-colon after an end statement.)

Checking for a unit first coefficient is bad, and silly. Testing for exact equality of floating point numbers is never a good thing to do. Instead, simply divide the coefficients by that first coefficient, (at least after testing that it is non-zero.) This works ALWAYS.

Next, if you are going to generate an error and not return a result, then use the function error! If you prefer, rather than disp use the warning function instead.

[Isaac Mancero Mosquera]

23/May/2012: There was an error in the "k2" expression in the file I submitted yesterday.
Now is corrected.