I coded these equations. The equation in the picture for s is the third root. The sgn(function) do I need to add this now?
I was not looking at the correctness of the equations merely at a solution for the equations.
All roots to a cubic are complex as the reals are a subset of the complex. If you do not make the variables complex you quickly get NAN, the set 1,1,1,1, returns t as NAN using reals. We can call it a real say X = 1, but it is in real terms X = 1 +0i.
I took complex theory at ANU in 1976 with a brilliant guy. We learnt a lot of interesting stuff.
As I said pure math is always QA.
The program is just a simple start, it is not meant to be exhaustive, nor did I test it beyond one trivial example. It was just a bit of fun, as is most pure math in the end for the first 200 years.
I cannot load a zip file here and I do not have github, I will put the code on Fortran Forum so you may play with it.
Humans did not invent imaginary numbers, they found them yes, but they did not invent them.
The program is simply wrong. Not one of the “roots” that it prints satisfies the cubic equation – not one of the output values x_1, x_2, x_3 satisfies ax^3+bx^2+cx+d = 0 for the values a=b=c=d=1 that you chose.
If I understand correctly, this is from your earlier post, please correct me if I am wrong. As it is stated that It has been verified I did not check it.
If I take all one’s as noted then Q = 3-1/9 which is 0.2222 as a real. So that is correct.
R is (9-27-2)/54 which is -0.3704, which is a real as found.
The next equation is sqrt((qqq) + (r*r)) = sqrt(0.01 + 0.1371) = 0.384
for S the absolute values of sqrt((qqq) + (r*r)) > R so S is positive real (1/3) root of 0.0144 which according to EXCEL is about 0.244 as shown in the program.
The next T is the real (1/3) root of (-0.754), which I coded very badly assuming that ** could handle negative numbers and I am wrong.
The required code is
if((R%RE-SD%RE) .lt. 0.0D0) then
T = -((abs(R-SD))**THIRDROOT)
else
T = (R - SD)**THIRDROOT
end if
The solution is a set of complex numbers (1,0),(0,1),(0,-1)
This is why we have QA. In a better version we would keep track of all issues like that, but ** why not overload **.
There are more checks to make, I was merely having some fun.
If you look at the first post in this thread, you will note that the topic concerns cubic equations with purely real coefficients. The earlier program listings contain mostly real variables, with just a few complex variables added as needed for the case where only one root is real.
The “verified” solution has the same context. It is the generalization to complex variables, coupled with using **(1/3.0) with a complex argument, that gave rise to the errors. See, for example, this paper, which stresses on the first page that taking the principal value is incorrect in the context of solving cubic equations.
I understand your points. I do not disagree with your points, but in Fortran you need to use Complex to solve even 1,1,1,1 problem as T is negative. I was merely looking as you are to solve real a,b,c,d - but you can set them up as complex, it is simple, so why not. I thought the answer looked strange, but I had a 15 yr old daughter making sounds that we were leaving so I just posted it. I knew you would catch the errors. Sorry, but you are great at that.
I did not put in limit tests, this site is for experts not for someone who cannot follow the math, well partly.
I had a few minutes to spare and I thought coding it would be fun, I am really not interested in the underlying problems, anyone who wants to can extend the program, I made it freely available.
It was just a bit of fun and any day where I converse with you is a good day.
I added a root checker to verify the results, as long as all the a,b,c,d > 0 or < 0 then the original equations work, but if you change one of the a,b,c,d to a negative only then the original equations do not work, not even close.
Plus you do not need a, it is overspecified, simply divide by a then b = x1+x2+x3 c = x1x2 + x1x3+x2x3
and d = x1x2x3 allowing for the sign.
