!*********************************************************** !* Solve complex equation of degree 3: * !* a z^3 + b z^2 + c z + d = 0 * !* ------------------------------------------------------- * !* SAMPLE RUN: * !* a = (1,0) b = (2,-10) c = (-4,1) d = (5,2) * !* * !* Equation az**3 + bz**2 + c*z + d = 0 * !* Solutions z: * !* z = (-0.4702985,0.6684887) * !* F(z) = (3.4763884E-06,-6.9430839E-06) * !* z = (0.6017774,-0.2911127) * !* F(z) = (1.1284958E-05,-2.6326748E-06) * !* z = (-2.131480,9.622624) * !* F(z) = (5.7570851E-05,5.4173665E-06) * !* * !* J-P Moreau, Paris. * !* (www.jpmoreau.fr) * !*********************************************************** ! NOTE: ! z = Z - (b/3a) ! p = (c/a) - (b²/3a²) ! q = (2b^3/27a^3) + (d/a) - (bc)/3a²)! ! ! Now the equation is reduced to Z^3 + pZ + q = 0 ! ! with the roots: ! ! Z[1..3] = cubroot(-q/2 + sqrt(q²/4 + p^3/27)) + cubroot(-q/2 - sqrt(q²/4 + p^3/27)) !------------------------------------------------------------------------------------ Program CRoot3 Complex a,b,c,d, det, p, q, sdet, u, u1, u2, v, v1, v2, w, z1 Complex z(9) ! seek three complex roots of complex equation: ! a z**3 + b z**2 + c*z + d = 0 a = CMPLX(1.0,0.0) b = CMPLX(2.0,-10.0) c = CMPLX(-4.0,1.0) d = CMPLX(5.0,2.0) ! calculate p = c/a - b²/3a² p = (c/a) - (b*b/(3*a*a)) ! calculate q = 2b**3/27a**3 + d/a - bc/3a² q = (2*b*b*b/(27*a*a*a)) + (d/a) - (b*c/(3*a*a)) ! calculate det = q²/4 + p^3/27 det = (q*q/4.) + (p*p*p/27.) sdet = SQRT(det) !now sdet contains sqrt(q²/4 + p**3/27 v = -q/2. + sdet call CUBRT(v,u,u1,u2) !3 cubic roots ! print *, u, u1, u2 w = -q/2. - sdet ! print *,' w=', w call CUBRT(w,v,v1,v2) !3 cubic roots ! print *, v, v1, v2 z(1) = u + v z(2) = u + v1 z(3) = u + v2 z(4) = u1 + v z(5) = u1 + v1 z(6) = u1 + v2 z(7) = u2 + v z(8) = u2 + v1 z(9) = u2 + v2 ! Note: only 3 z(i) are correct solutions. do i=1, 9 z(i) = z(i) - (b/(3*a)) end do print *,' ' print *,' Equation az**3 + bz**2 + c*z + d = 0' print *,' Solutions z:' do i=1, 9 call F(a,b,c,d,z(i),z1) if (ABS(z1) < 0.001) then !detect if solution is correct write(*,*) ' z = ', z(i); write(*,*) ' F(z) = ', z1 end if end do print *,' ' stop END Function ATAN1(xnumerator, denominator) ! Return a phase between -PI and +PI PI = 4.0*atan(1.0) if (abs(denominator) < 1E-15) then if (xnumerator < 0.0) then ATAN1 = -PI/2.0 else ATAN1 = PI/2.0 end if else phase=atan(xnumerator/denominator) if (denominator < 0.0) phase = phase + PI ATAN1 = phase end if return End !ATAN1 ! z1, z2, z3 = z**(1/3) Subroutine CUBRT(z, z1, z2, z3) Complex z,z1,z2,z3 PI = 4.0*atan(1.0) r=abs(z) r=r**(1./3.) if (ABS(REAL(z))<1.E-15) then if (AIMAG(z) < 0.) then t = -PI/2. else t = PI/2. end if else ! return t between -pi and +pi t=ATAN1(AIMAG(z),REAL(z)) end if t=t/3. z1=CMPLX(r*cos(t),r*sin(t)) tt=t-(2*PI/3.) z2=CMPLX(r*cos(tt),r*sin(tt)) tt=t+(2*PI/3); z3=CMPLX(r*cos(tt),r*sin(tt)) return End ! z1 = a*z^3 + b*z2 + c*z + d Subroutine F(a,b,c,d,z,z1) Complex a,b,c,d,z,z1 z1 = a*z*z*z + b*z*z + c*z + d return End !end of file croot3.f90