Skip to content

Commit

Permalink
Adding complex linear/quadratic/cubic equation roots
Browse files Browse the repository at this point in the history
  • Loading branch information
@jonathanolson
jonathanolson committed Jul 14, 2023
1 parent 61402ec commit ff2cf28
Show file tree
Hide file tree
Showing 2 changed files with 261 additions and 2 deletions.
162 changes: 162 additions & 0 deletions js/Complex.js
View file
Original file line number Diff line number Diff line change
Expand Up @@ -85,6 +85,20 @@ class Complex {
return this.getMagnitudeSquared();
}

/**
* Returns the argument of this complex number (immutable)
* @public
*
* @returns {number}
*/
getArgument() {
return Math.atan2( this.imaginary, this.real );
}

get argument() {
return this.getArgument();
}

/**
* Exact equality comparison between this Complex and another Complex.
* @public
Expand Down Expand Up @@ -168,6 +182,17 @@ class Complex {
);
}

/**
* Complex negation
* Immutable version of negate
* @public
*
* @returns {Complex}
*/
negated() {
return new Complex( -this.real, -this.imaginary );
}

/**
* Square root.
* Immutable form of sqrt.
Expand Down Expand Up @@ -380,6 +405,16 @@ class Complex {
);
}

/**
* Mutable Complex negation
* @public
*
* @returns {Complex}
*/
negate() {
return this.setRealImaginary( -this.real, -this.imaginary );
}

/**
* Sets this Complex to e to the power of this complex number. $e^{a+bi}=e^a\cos b + i\sin b$.
* This is the mutable version of exponentiated
Expand Down Expand Up @@ -455,6 +490,27 @@ class Complex {
return this.setRealImaginary( this.real, -this.imaginary );
}

/**
* Returns the cube roots of this complex number.
* @public
*
* @returns {Complex[]}
*/
getCubeRoots() {
const arg3 = this.argument / 3;
const abs = this.magnitude;

const really = Complex.real( Math.cbrt( abs ) );

const principal = really.times( Complex.imaginary( arg3 ).exponentiate() );

return [
principal,
really.times( Complex.imaginary( arg3 + Math.PI * 2 / 3 ).exponentiate() ),
really.times( Complex.imaginary( arg3 - Math.PI * 2 / 3 ).exponentiate() )
];
}

/**
* Debugging string for the complex number (provides real and imaginary parts).
* @public
Expand Down Expand Up @@ -499,6 +555,112 @@ class Complex {
static createPolar( magnitude, phase ) {
return new Complex( magnitude * Math.cos( phase ), magnitude * Math.sin( phase ) );
}

/**
* Returns an array of the roots of the quadratic equation $ax + b=0$, or null if every value is a solution.
* @public
*
* @param {Complex} a
* @param {Complex} b
* @returns {Array.<Complex>|null} - The roots of the equation, or null if all values are roots.
*/
static solveLinearRoots( a, b ) {
if ( a.equals( Complex.ZERO ) ) {
return b.equals( Complex.ZERO ) ? null : [];
}

return [ b.dividedBy( a ).negate() ];
}

/**
* Returns an array of the roots of the quadratic equation $ax^2 + bx + c=0$, or null if every value is a
* solution.
* @public
*
* @param {Complex} a
* @param {Complex} b
* @param {Complex} c
* @returns {Array.<Complex>|null} - The roots of the equation, or null if all values are roots (if multiplicity>1,
* returns multiple copies)
*/
static solveQuadraticRoots( a, b, c ) {
if ( a.equals( Complex.ZERO ) ) {
return Complex.solveLinearRoots( b, c );
}

const denom = Complex.real( 2 ).multiply( a );
const d1 = b.times( b );
const d2 = Complex.real( 4 ).multiply( a ).multiply( c );
const discriminant = d1.subtract( d2 ).sqrt();
return [
discriminant.minus( b ).divide( denom ),
discriminant.negated().subtract( b ).divide( denom )
];
}

/**
* Returns an array of the roots of the cubic equation $ax^3 + bx^2 + cx + d=0$, or null if every value is a
* solution.
* @public
*
* @param {Complex} a
* @param {Complex} b
* @param {Complex} c
* @param {Complex} d
* @returns {Array.<Complex>|null} - The roots of the equation, or null if all values are roots (if multiplicity>1,
* returns multiple copies)
*/
static solveCubicRoots( a, b, c, d ) {
if ( a.equals( Complex.ZERO ) ) {
return Complex.solveQuadraticRoots( b, c, d );
}

const denom = a.times( Complex.real( 3 ) ).negate();
const a2 = a.times( a );
const b2 = b.times( b );
const b3 = b2.times( b );
const c2 = c.times( c );
const c3 = c2.times( c );
const abc = a.times( b ).times( c );

// TODO: factor out constant numeric values

const D0_1 = b2;
const D0_2 = a.times( c ).times( Complex.real( 3 ) );
const D1_1 = b3.times( Complex.real( 2 ) ).add( a2.times( d ).multiply( Complex.real( 27 ) ) );
const D1_2 = abc.times( Complex.real( 9 ) );

if ( D0_1.equals( D0_2 ) && D1_1.equals( D1_2 ) ) {
const tripleRoot = b.divide( denom );
return [ tripleRoot, tripleRoot, tripleRoot ];
}

const Delta0 = D0_1.minus( D0_2 );
const Delta1 = D1_1.minus( D1_2 );

const discriminant1 = abc.times( d ).multiply( Complex.real( 18 ) ).add( b2.times( c2 ) );
const discriminant2 = b3.times( d ).multiply( Complex.real( 4 ) )
.add( c3.times( a ).multiply( Complex.real( 4 ) ) )
.add( a2.times( d ).multiply( d ).multiply( Complex.real( 27 ) ) );

if ( discriminant1.equals( discriminant2 ) ) {
const simpleRoot = (
abc.times( Complex.real( 4 ) ).subtract( b3.plus( a2.times( d ).multiply( Complex.real( 9 ) ) ) )
).divide( a.times( Delta0 ) );
const doubleRoot = ( a.times( d ).multiply( Complex.real( 9 ) ).subtract( b.times( c ) ) ).divide( Delta0.times( Complex.real( 2 ) ) );
return [ simpleRoot, doubleRoot, doubleRoot ];
}
let Ccubed;
if ( D0_1.equals( D0_2 ) ) {
Ccubed = Delta1;
}
else {
Ccubed = Delta1.plus( ( Delta1.times( Delta1 ).subtract( Delta0.times( Delta0 ).multiply( Delta0 ).multiply( Complex.real( 4 ) ) ) ).sqrt() ).divide( Complex.real( 2 ) );
}
return Ccubed.getCubeRoots().map( root => {
return b.plus( root ).add( Delta0.dividedBy( root ) ).divide( denom );
} );
}
}

dot.register( 'Complex', Complex );
Expand Down
101 changes: 99 additions & 2 deletions js/ComplexTests.js
View file
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
// Copyright 2017-2022, University of Colorado Boulder

/**
* Bounds2 tests
* Complex.js tests
*
* @author Jonathan Olson (PhET Interactive Simulations)
* @author Sam Reid (PhET Interactive Simulations)
Expand Down Expand Up @@ -63,4 +63,101 @@ QUnit.test( 'Sin of', assert => {
const a = new Complex( 1, 1 );
const b = new Complex( 1.29845758, 0.634963914 );
approximateComplexEquals( assert, a.sinOf(), b, 'Sin Of' );
} );
} );

QUnit.test( 'getCubeRoots', assert => {
const actual = new Complex( 0, 8 ).getCubeRoots();
const expected = [
new Complex( Math.sqrt( 3 ), 1 ),
new Complex( -Math.sqrt( 3 ), 1 ),
new Complex( 0, -2 )
];
approximateComplexEquals( assert, expected[ 0 ], actual[ 0 ], 'root 1' );
approximateComplexEquals( assert, expected[ 1 ], actual[ 1 ], 'root 2' );
approximateComplexEquals( assert, expected[ 2 ], actual[ 2 ], 'root 3' );
} );

QUnit.test( 'linear roots', assert => {
approximateComplexEquals( assert, Complex.real( -2 ), Complex.solveLinearRoots( Complex.real( 3 ), Complex.real( 6 ) )[ 0 ], '3x + 6 = 0 // x=-2' );
approximateComplexEquals( assert, new Complex( -2, -1 ), Complex.solveLinearRoots( Complex.real( 3 ), new Complex( 6, 3 ) )[ 0 ], '3x + 6 + 3i = 0 // x=-2-i' );
approximateComplexEquals( assert, Complex.real( -3 ), Complex.solveLinearRoots( new Complex( 2, 1 ), new Complex( 6, 3 ) )[ 0 ], '(2 + i)x + 6 + 3i = 0 // x=-3' );
} );

QUnit.test( 'quadratic roots 0', assert => {
const roots0 = Complex.solveQuadraticRoots( Complex.real( 1 ), Complex.real( -3 ), Complex.real( 2 ) );
const roots0_0 = Complex.real( 2 );
const roots0_1 = Complex.real( 1 );
approximateComplexEquals( assert, roots0[ 0 ], roots0_0, 'x^2 - 3x + 2 = 0 // x=2 case' );
approximateComplexEquals( assert, roots0[ 1 ], roots0_1, 'x^2 - 3x + 2 = 0 // x=1 case' );
} );

QUnit.test( 'quadratic roots 1', assert => {
const roots1 = Complex.solveQuadraticRoots( Complex.real( 2 ), Complex.real( 8 ), Complex.real( 8 ) );
const roots1_x = Complex.real( -2 );
approximateComplexEquals( assert, roots1[ 0 ], roots1_x, '8x^2 + 8x + 2 = 0 // x=-2 case (first)' );
approximateComplexEquals( assert, roots1[ 1 ], roots1_x, '8x^2 + 8x + 2 = 0 // x=-2 case (second)' );
} );

QUnit.test( 'quadratic roots 2', assert => {
const roots2 = Complex.solveQuadraticRoots( Complex.real( 1 ), Complex.real( 0 ), Complex.real( -2 ) );
const roots2_0 = Complex.real( Math.sqrt( 2 ) );
const roots2_1 = Complex.real( -Math.sqrt( 2 ) );
approximateComplexEquals( assert, roots2[ 0 ], roots2_0, 'x^2 - 2 = 0 // x=sqrt(2) case' );
approximateComplexEquals( assert, roots2[ 1 ], roots2_1, 'x^2 - 2 = 0 // x=-sqrt(2) case' );
} );

QUnit.test( 'quadratic roots 3', assert => {
const roots3 = Complex.solveQuadraticRoots( Complex.real( 1 ), Complex.real( -2 ), Complex.real( 2 ) );
const roots3_0 = new Complex( 1, 1 );
const roots3_1 = new Complex( 1, -1 );
approximateComplexEquals( assert, roots3[ 0 ], roots3_0, 'x^2 - 2x + 2 = 0 // x=1+i case' );
approximateComplexEquals( assert, roots3[ 1 ], roots3_1, 'x^2 - 2x + 2 = 0 // x=1-i case' );
} );

QUnit.test( 'quadratic roots 4', assert => {
const roots4 = Complex.solveQuadraticRoots( Complex.real( 1 ), new Complex( -3, -2 ), new Complex( 1, 3 ) );
const roots4_0 = new Complex( 2, 1 );
const roots4_1 = new Complex( 1, 1 );
approximateComplexEquals( assert, roots4[ 0 ], roots4_0, '(1 + 3i)x^2 + (-3 - 2i)x + 1 = 0 // x=2+i case' );
approximateComplexEquals( assert, roots4[ 1 ], roots4_1, '(1 + 3i)x^2 + (-3 - 2i)x + 1 = 0 // x=1+i case' );
} );

QUnit.test( 'cubic roots 0', assert => {
const roots = Complex.solveCubicRoots( Complex.real( 1 ), Complex.real( -6 ), Complex.real( 11 ), Complex.real( -6 ) );
const roots_0 = Complex.real( 1 );
const roots_1 = Complex.real( 3 );
const roots_2 = Complex.real( 2 );
approximateComplexEquals( assert, roots[ 0 ], roots_0, 'x^3 - 6x^2 + 11x - 6 = 0 // x=1 case' );
approximateComplexEquals( assert, roots[ 1 ], roots_1, 'x^3 - 6x^2 + 11x - 6 = 0 // x=3 case' );
approximateComplexEquals( assert, roots[ 2 ], roots_2, 'x^3 - 6x^2 + 11x - 6 = 0 // x=2 case' );
} );

QUnit.test( 'cubic roots 1', assert => {
const roots = Complex.solveCubicRoots( Complex.real( 1 ), Complex.real( 0 ), Complex.real( 0 ), Complex.real( -8 ) );
const roots_0 = new Complex( -1, -Math.sqrt( 3 ) );
const roots_1 = Complex.real( 2 );
const roots_2 = new Complex( -1, Math.sqrt( 3 ) );
approximateComplexEquals( assert, roots[ 0 ], roots_0, 'x^3 - 8 = 0 // x=-1-sqrt(3)i case' );
approximateComplexEquals( assert, roots[ 1 ], roots_1, 'x^3 - 8 = 0 // x=2 case' );
approximateComplexEquals( assert, roots[ 2 ], roots_2, 'x^3 - 8 = 0 // x=-1+sqrt(3)i case' );
} );

QUnit.test( 'cubic roots 2', assert => {
const roots = Complex.solveCubicRoots( Complex.real( 2 ), Complex.real( 8 ), Complex.real( 8 ), Complex.real( 0 ) );
const roots_0 = Complex.real( 0 );
const roots_1 = Complex.real( -2 );
const roots_2 = Complex.real( -2 );
approximateComplexEquals( assert, roots[ 0 ], roots_0, '2x^3 + 8x^2 + 8x = 0 // x=0' );
approximateComplexEquals( assert, roots[ 1 ], roots_1, '2x^3 + 8x^2 + 8x = 0 // x=-2 case' );
approximateComplexEquals( assert, roots[ 2 ], roots_2, '2x^3 + 8x^2 + 8x = 0 // x=-2 case' );
} );

QUnit.test( 'cubic roots 3', assert => {
const roots = Complex.solveCubicRoots( Complex.real( 1 ), Complex.real( 1 ), Complex.real( 1 ), Complex.real( 1 ) );
const roots_0 = Complex.real( -1 );
const roots_1 = new Complex( 0, -1 );
const roots_2 = new Complex( 0, 1 );
approximateComplexEquals( assert, roots[ 0 ], roots_0, 'x^3 + x^2 + x + 1 = 0 // x=-1' );
approximateComplexEquals( assert, roots[ 1 ], roots_1, 'x^3 + x^2 + x + 1 = 0 // x=-i case' );
approximateComplexEquals( assert, roots[ 2 ], roots_2, 'x^3 + x^2 + x + 1 = 0 // x=i case' );
} );

0 comments on commit ff2cf28

Please sign in to comment.