Solving cubic polynomials - Revision history

Jpye: /* Calculation checking */

2022-07-07T02:50:07Z

Calculation checking

← Older revision		Revision as of 02:50, 7 July 2022
Line 183:		Line 183:
	The following [https://docs.sympy.org/latest/tutorial/simplification.html sympy] code can be used to confirm the algebra of the 'general solution' above. Note that this page uses the depressed form <math>y^3-3Py-2Q=0</math>, while most other sources use <math>y^3+py +q=0</math>, with results that are equivalent but have more constants lying around.		The following [https://docs.sympy.org/latest/tutorial/simplification.html sympy] code can be used to confirm the algebra of the 'general solution' above. Note that this page uses the depressed form <math>y^3-3Py-2Q=0</math>, while most other sources use <math>y^3+py +q=0</math>, with results that are equivalent but have more constants lying around.

	To determine the expressions for <math>P</math> and <math>Q</math> in the depressed cubic. This ~~expression with~~ <math>~~a_3~~ = a/3</math> ~~aims~~ at ~~improved~~ numerical precision.		To determine the expressions for <math>P</math> and <math>Q</math> in the depressed cubic. This unnecessary replacement <code>a3</code><math>= a/3</math> is aimed at increasing numerical precision.
	<source lang="py">		<source lang="py">
	from sympy import *		from sympy import *

Jpye: /* Implementation in C */

2022-07-07T02:47:36Z

Implementation in C

← Older revision		Revision as of 02:47, 7 July 2022
Line 168:		Line 168:

	<source lang=py>		<source lang=py>
	from sympy import *; x,a,b,c,F = symbols('x,a,b,c,F')		from sympy import *; x,a,b,c,w,F = symbols('x,a,b,c,w,F')
	cubic = x*3 + ax*2 + bx + c		cubic = x*3 + ax*2 + bx + c
	(cubic.subs(x,w*F)/F^3).expand()		(cubic.subs(x,w*F)/F^3).expand()
	</source>		</source>

Jpye: /* Implementation in C */

2022-07-07T02:46:38Z

Implementation in C

← Older revision		Revision as of 02:46, 7 July 2022
Line 154:		Line 154:
	This approach has been implemented in plain C code as {{srcbranch\|python3\|models/johnpye/fprops/cubicroots.c}}, licensed under the GPL. The code includes debugging output that be turned off by disabling <code>CUBICROOTS_DEBUG</code>. The routine can be compiled to calculate using standard <code>double</code> precision, or with the switches <code>CUBICROOTS_EXTRA_PRECISION</code> and <code>CUBICROOTS_LOW_PRECISION</code>, calculation with <code>long double</code> or <code>float</code>, respectively, can be selected. This selection is done at compile time only, and only affects the internal precision of the calculation, not the API. Currently this code is embedded within [[FPROPS]] build (which uses [https://scons.org SCons]), but it is trivial to convert it to a standalone code. The dependencies on <code>color.h</code> and <code>common.h</code> can very easily be removed.		This approach has been implemented in plain C code as {{srcbranch\|python3\|models/johnpye/fprops/cubicroots.c}}, licensed under the GPL. The code includes debugging output that be turned off by disabling <code>CUBICROOTS_DEBUG</code>. The routine can be compiled to calculate using standard <code>double</code> precision, or with the switches <code>CUBICROOTS_EXTRA_PRECISION</code> and <code>CUBICROOTS_LOW_PRECISION</code>, calculation with <code>long double</code> or <code>float</code>, respectively, can be selected. This selection is done at compile time only, and only affects the internal precision of the calculation, not the API. Currently this code is embedded within [[FPROPS]] build (which uses [https://scons.org SCons]), but it is trivial to convert it to a standalone code. The dependencies on <code>color.h</code> and <code>common.h</code> can very easily be removed.

	There is also a test suite that can be compiled by passing <code>-DTEST</code> to the compiler. The various test cases are implemented via the <code>RTEST(a,b,c,nr,x0,x1,x2)</code> macro, which checks the that the roots of the polynomial <math>x^3 + a x^2 + b x + c = 0</math> solve to the expected <code>nr</code> roots, 1 or 3, with values <math>\left\lbrace x_0\right\rbrace</math> or <math>\left\lbrace x_0, x_1, x_2\right\rbrace</math>, respectively. The solved roots are checked to be correct with relative tolerance of <code>TESTTOL</code>, which is nominally <code>1-e13</code> (at <~~math~~>double</~~math~~> precision). Further testing may reveal that this tolerance is overly optimistic, but it works so far.		There is also a test suite that can be compiled by passing <code>-DTEST</code> to the compiler. The various test cases are implemented via the <code>RTEST(a,b,c,nr,x0,x1,x2)</code> macro, which checks the that the roots of the polynomial <math>x^3 + a x^2 + b x + c = 0</math> solve to the expected <code>nr</code> roots, 1 or 3, with values <math>\left\lbrace x_0\right\rbrace</math> or <math>\left\lbrace x_0, x_1, x_2\right\rbrace</math>, respectively. The solved roots are checked to be correct with relative tolerance of <code>TESTTOL</code>, which is nominally <code>1-e13</code> (at <code>double</code> precision). Further testing may reveal that this tolerance is overly optimistic, but it works so far.

	Some challenging cases were found during testing, such as		Some challenging cases were found during testing, such as

Jpye at 02:45, 7 July 2022

2022-07-07T02:45:39Z

← Older revision		Revision as of 02:45, 7 July 2022
Line 152:		Line 152:
	== Implementation in C ==		== Implementation in C ==

	This approach has been implemented in plain C code as {{srcbranch\|python3\|models/johnpye/fprops/cubicroots.c}}, licensed under the GPL. The code includes debugging output that be turned off by disabling <code>CUBICROOTS_DEBUG</code>. The routine can be compiled to calculate using standard <code>double</code> precision, or with the switches <code>CUBICROOTS_EXTRA_PRECISION</code> and <code>CUBICROOTS_LOW_PRECISION</code>, calculation with <code>long double</code> or <code>float</code>, respectively, can be selected. This selection is done at compile time only, and only affects the internal precision of the calculation, not the API. Currently this code is embedded within [[FPROPS]] build (which uses [https://scons.org] SCons]), but it is trivial to convert it to a standalone code. The dependencies on <code>color.h</code> and <code>common.h</code> can very easily be removed.		This approach has been implemented in plain C code as {{srcbranch\|python3\|models/johnpye/fprops/cubicroots.c}}, licensed under the GPL. The code includes debugging output that be turned off by disabling <code>CUBICROOTS_DEBUG</code>. The routine can be compiled to calculate using standard <code>double</code> precision, or with the switches <code>CUBICROOTS_EXTRA_PRECISION</code> and <code>CUBICROOTS_LOW_PRECISION</code>, calculation with <code>long double</code> or <code>float</code>, respectively, can be selected. This selection is done at compile time only, and only affects the internal precision of the calculation, not the API. Currently this code is embedded within [[FPROPS]] build (which uses [https://scons.org SCons]), but it is trivial to convert it to a standalone code. The dependencies on <code>color.h</code> and <code>common.h</code> can very easily be removed.

	There is also a test suite that can be compiled by passing <code>-DTEST</code> to the compiler. The various test cases are implemented via the <code>RTEST(a,b,c,nr,x0,x1,x2)</code> macro, which checks the that the roots of the polynomial <math>x^3 + a x^2 + b x + c = 0</math> solve to the expected <code>nr</code> roots, 1 or 3, with values <math>\left\lbrace x_0\right\rbrace</math> or <math>\left\lbrace x_0, x_1, x_2\right\rbrace</math>, respectively. The solved roots are checked to be correct with relative tolerance of <code>TESTTOL</code>, which is nominally <code>1-e13</code> (at <math>double</math> precision). Further testing may reveal that this tolerance is overly optimistic, but it works so far.		There is also a test suite that can be compiled by passing <code>-DTEST</code> to the compiler. The various test cases are implemented via the <code>RTEST(a,b,c,nr,x0,x1,x2)</code> macro, which checks the that the roots of the polynomial <math>x^3 + a x^2 + b x + c = 0</math> solve to the expected <code>nr</code> roots, 1 or 3, with values <math>\left\lbrace x_0\right\rbrace</math> or <math>\left\lbrace x_0, x_1, x_2\right\rbrace</math>, respectively. The solved roots are checked to be correct with relative tolerance of <code>TESTTOL</code>, which is nominally <code>1-e13</code> (at <math>double</math> precision). Further testing may reveal that this tolerance is overly optimistic, but it works so far.

Jpye at 02:38, 7 July 2022

2022-07-07T02:38:12Z

← Older revision		Revision as of 02:38, 7 July 2022
Line 1:		Line 1:
	The purpose of this page is to describe a robust process, in software, for finding real roots of a real cubic polynomial, something that is essential for dealing with cubic equations of state in chemical process modelling. It turns out that this is one of the great problems in the history of mathematics, and not that simple, as [https://youtu.be/N-KXStupwsc this video] helpfully explains. This page attempts to stick to the essential details and describe how to implement a solution suitable for coding.		The purpose of this page is to describe a robust process, in software, for finding real roots of a real cubic polynomial, something that is essential for dealing with cubic equations of state in chemical process modelling. It turns out that this is one of the great problems in the history of mathematics, and not that simple, as [https://youtu.be/N-KXStupwsc this video from Burkard Polster] helpfully explains. This page attempts to stick to the essential details and describe how to implement a solution suitable for coding.

	== General solution ==		== General solution ==

Jpye at 02:37, 7 July 2022

2022-07-07T02:37:55Z

← Older revision		Revision as of 02:37, 7 July 2022
Line 1:		Line 1:
	The purpose of this page is to describe a robust process, in software, for finding real roots of a real cubic polynomial, something that is essential for dealing with cubic equations of state in chemical process modelling. It turns out that this is one of the great problems in the history of mathematics, and not ~~at all~~ simple, as [https://youtu.be/N-KXStupwsc this video] helpfully explains. This page attempts to stick to the essential details and describe how to implement a solution suitable for coding.		The purpose of this page is to describe a robust process, in software, for finding real roots of a real cubic polynomial, something that is essential for dealing with cubic equations of state in chemical process modelling. It turns out that this is one of the great problems in the history of mathematics, and not that simple, as [https://youtu.be/N-KXStupwsc this video] helpfully explains. This page attempts to stick to the essential details and describe how to implement a solution suitable for coding.

	== General solution ==		== General solution ==

Jpye at 02:37, 7 July 2022

2022-07-07T02:37:37Z

← Older revision		Revision as of 02:37, 7 July 2022
Line 1:		Line 1:
	The purpose of this page is to describe a robust process, in software, for finding real roots of a real cubic polynomial. It turns out that this is one of the great problems in the history of mathematics, and not at all simple, as [https://youtu.be/N-KXStupwsc this video] helpfully explains. This page attempts to stick to the essential details and describe how to implement a solution suitable for coding.		The purpose of this page is to describe a robust process, in software, for finding real roots of a real cubic polynomial, something that is essential for dealing with cubic equations of state in chemical process modelling. It turns out that this is one of the great problems in the history of mathematics, and not at all simple, as [https://youtu.be/N-KXStupwsc this video] helpfully explains. This page attempts to stick to the essential details and describe how to implement a solution suitable for coding.

	== General solution ==		== General solution ==

Jpye: /* Implementation in C */

2022-07-06T11:27:12Z

Implementation in C

← Older revision		Revision as of 11:27, 6 July 2022
Line 173:		Line 173:
	</source>		</source>

	Some subtle logic was also needed for cases where <math>P</math> and <math>D</math> are close to zero. The variable <code>TOLD</code> is used to check how small <math>\left\|D\right\|</math> can be before it is assumed to be zero, yielding double and single real roots, as opposed to a single or three different real roots. It is not clear yet whether this value of <code>TOLD</code> has been wisely chosen; for now it works with the test suite provided. The tolerances ~~<code>TOLQ</code> (~~for <math>Q</math>~~) and~~ <code>~~TOLP~~</code> (for <math>P</math>) work similarly.		Some subtle logic was also needed for cases where <math>P</math> and <math>D</math> are close to zero. The variable <code>TOLD</code> is used to check how small <math>\left\|D\right\|</math> can be before it is assumed to be zero, yielding double and single real roots, as opposed to a single or three different real roots. It is not clear yet whether this value of <code>TOLD</code> has been wisely chosen; for now it works with the test suite provided, which includes all of the test cases used in GSL as well as quite a few more. The tolerances for <math>Q</math> (<code>TOLQ</code>) and for <math>P</math> (<code>TOLP</code>) work similarly. Again, the values have not yet been carefully tuned.

	The logic for specifically handling the cases of <math>P=0</math> and <math>Q=0</math> adds a bit of code, but would reduce computation cost in cases where a lot of such polynomials are being solved. Perhaps it could be removed to make the code simpler.		The logic for specifically handling the cases of <math>P=0</math> and <math>Q=0</math> adds a bit of code, but would presumably reduce computation cost in cases where a lot of such polynomials are being solved. Perhaps it could be removed to make the code simpler.

	A current limitation of the code is that it has no way of passing back the nature of the roots, for example if the value of <math>D</math> is close to zero, or whether there were multiple roots or not. Only the number of roots is returned. It is not yet know if users of this function could benefit from that information or not. Note that the root are sorted so as to be returned in increasing order. Double roots could be the first two or the second two roots.		A current limitation of the code is that it has no way of passing back the nature of the roots, for example if the value of <math>D</math> is close to zero, or whether there were multiple roots or not. Only the number of roots is returned. It is not yet know if users of this function could benefit from that information or not. Note that the root are sorted so as to be returned in increasing order. Double roots could be the first two or the second two roots.

Jpye: /* Implementation in C */

2022-07-05T04:37:13Z

Implementation in C

← Older revision		Revision as of 04:37, 5 July 2022
Line 166:		Line 166:

	Both of these cases were improved by adding some additional logic for re-scaling <math>w = \frac{x}{F}</math>, and then also dividing the parameters <math>a,b,c</math> in order to keep the polynomial monic.		Both of these cases were improved by adding some additional logic for re-scaling <math>w = \frac{x}{F}</math>, and then also dividing the parameters <math>a,b,c</math> in order to keep the polynomial monic.

			<source lang=py>
			from sympy import *; x,a,b,c,F = symbols('x,a,b,c,F')
			cubic = x*3 + ax*2 + bx + c
			(cubic.subs(x,w*F)/F^3).expand()
			</source>

	Some subtle logic was also needed for cases where <math>P</math> and <math>D</math> are close to zero. The variable <code>TOLD</code> is used to check how small <math>\left\|D\right\|</math> can be before it is assumed to be zero, yielding double and single real roots, as opposed to a single or three different real roots. It is not clear yet whether this value of <code>TOLD</code> has been wisely chosen; for now it works with the test suite provided. The tolerances <code>TOLQ</code> (for <math>Q</math>) and <code>TOLP</code> (for <math>P</math>) work similarly.		Some subtle logic was also needed for cases where <math>P</math> and <math>D</math> are close to zero. The variable <code>TOLD</code> is used to check how small <math>\left\|D\right\|</math> can be before it is assumed to be zero, yielding double and single real roots, as opposed to a single or three different real roots. It is not clear yet whether this value of <code>TOLD</code> has been wisely chosen; for now it works with the test suite provided. The tolerances <code>TOLQ</code> (for <math>Q</math>) and <code>TOLP</code> (for <math>P</math>) work similarly.

Jpye: /* Implementation in C */

2022-07-03T11:53:24Z

Implementation in C

← Older revision		Revision as of 11:53, 3 July 2022
Line 171:		Line 171:
	The logic for specifically handling the cases of <math>P=0</math> and <math>Q=0</math> adds a bit of code, but would reduce computation cost in cases where a lot of such polynomials are being solved. Perhaps it could be removed to make the code simpler.		The logic for specifically handling the cases of <math>P=0</math> and <math>Q=0</math> adds a bit of code, but would reduce computation cost in cases where a lot of such polynomials are being solved. Perhaps it could be removed to make the code simpler.

	A current limitation of the code is that it has no way of passing back the nature of the roots, for example if the value of <math>D</math> is close to zero, or whether there were multiple roots or not. Only the number of roots is returned. It is not yet know if users of this function could benefit from that information or not.		A current limitation of the code is that it has no way of passing back the nature of the roots, for example if the value of <math>D</math> is close to zero, or whether there were multiple roots or not. Only the number of roots is returned. It is not yet know if users of this function could benefit from that information or not. Note that the root are sorted so as to be returned in increasing order. Double roots could be the first two or the second two roots.

	== Calculation checking ==		== Calculation checking ==