!*************************************************************** !* Solve complex equation of degree 4: * !* a z^4 + b z^3 + c z^2 + d z + e = 0 * !* ----------------------------------------------------------- * !* SAMPLE RUN: * !* a = (1,0) b = (2,-10) c = (-4,1) d = (5,2) e = (3,-7.5) * !* * !* Equation az^4 + bz^3 + cz^2 + d z + e = 0 * !* Solutions z1..z4: * !* * !* z1 = (-2.139849,9.625143) * !* F(z1) = (-1.0131565E-04,1.4495029E-04) * !* z2 = (0.6901917,-0.7063792) * !* F(z2) = (1.9168254E-05,-2.2552895E-05) * !* z3 = (-0.8758767,0.2231503) * !* F(z3) = (2.3990218E-05,-1.5628273E-05) * !* z4 = (0.3255335,0.8580866) * !* F(z4) = (2.4645762E-05,-1.7702958E-05) * !* * !* ----------------------------------------------------------- * !* Reference: Mathématiques et statitiques - Programmes en * !* BASIC. Editions du P.S.I., 1981. * !* * !* Adapted to complex Domain By J-P Moreau, Paris. * !* (www.jpmoreau.fr) * !*************************************************************** ! Explanations: ! ------------ ! Let us solve complex equation z^4 + a*z^3 + b*z^2 + c*z + d = 0 (1) ! ! Let us set z = z1 - a/4 and equation (1) becomes: ! ! (z1²+2ckz1+l)(z1²-2ckz1+m) = 0 (2) ! ! where: l + m -4ck² = q | q = b - 3a²/8 ! 2 ck (m-1) = r | (3) with: r = c-(ab/2)+(a^3/8) ! lm = s | s = d-(ac/4)+(a²b/16)-(3a^4/256) ! ! ck² = zz is a solution of complex cubic equation: ! ! zz^3 + aa zz^2 + bb z + cc = 0 (4) ! ! with: aa = q/2, bb = (q²-4s)/16, cc = -r²/64 ! ! So we solve equation (4) using Croot3 (3 roots zz1, zz2, zz3), then ! we calculate complex l and m by solving system (3), finally, we solve ! both 2nd degree equations of product (2). We obtain three sets of four ! complex roots, only one of which is valid, i.e. annulates equation (1). ! !----------------------------------------------------------------------- Program Test_CRoot4 Complex a,b,c,d,e,z1,z2,z3,z4,zz print *,' ' print *,' Equation az^4 + bz^3 + cz^2 + d z + e = 0' print *,' Solutions z1..z4:' a = CMPLX(1.0, 0.0) b = CMPLX(2.0,-10.0) c = CMPLX(-4.0,1.0) d = CMPLX(5.0,2.0) e = CMPLX(3.0,-7.5) call Croot4(a,b,c,d,e, z1,z2,z3,z4) ! roots are in z1..z4, a,b,c,d are modified ! restore original a,b,c,d,e a = CMPLX(1.0, 0.0) b = CMPLX(2.0,-10.0) c = CMPLX(-4.0,1.0) d = CMPLX(5.0,2.0) e = CMPLX(3.0,-7.5) print *,' ' print *,' z1 = ', z1 call F4(a,b,c,d,e,z1,zz) print *,' F(z1) = ', zz print *,' z2 = ', z2 call F4(a,b,c,d,e,z2,zz) print *,' F(z2) = ', zz print *,' z3 = ', z3 call F4(a,b,c,d,e,z3,zz) print *,' F(z3) = ', zz print *,' z4 = ', z4 call F4(a,b,c,d,e,z4,zz) print *,' F(z4) = ', zz 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 F3(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 ! z1 = a*z^4+b*z^3+c*z^2+d^z+e Subroutine F4(a,b,c,d,e,z,z1) Complex a,b,c,d,e,z,z1 z1 = a*z**4 + b*z**3 + c*z*z + d*z + e return End ! z1 = z^3 + p*z + q (called by Croot3) !Subroutine F1(p, q, z, z1) ! Complex p,q,z,z1 ! z1 = z*z*z + p*z + q ! return !End ! z1 = a*z^3 + b*z2 + c*z + d Subroutine F1(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 ! z1 = z^4+a*z^3+b*z^2+c^z+d (called by Croot4) Subroutine F(a,b,c,d,z, z1) Complex a,b,c,d,z,z1 z1 = z**4 + a*z**3 + b*z*z + c*z + d return End ! Solve equation az*z+b*z+c = 0 ! return 2 solutions s1, s2. Subroutine Equa2c(a,b,c,s1,s2) COMPLEX a,b,c,s1,s2 REAL*8 r,sx,sy,ux,uy,vx,vy,wx,wy ux = REAL(b)*REAL(b) - AIMAG(b)*AIMAG(b) - 4.d0 * REAL(a)*REAL(c) + 4.d0 * AIMAG(a)*AIMAG(c) uy = 2.d0*REAL(b)*AIMAG(b) - 4.d0 * REAL(a)*AIMAG(c) - 4.d0*AIMAG(a)*REAL(c) r = dsqrt(ux*ux + uy*uy) vx = dsqrt((r+ux)/2.d0) vy = dsqrt((r-ux)/2.d0) if (uy<0.) vy = -vy wx = (-REAL(b) - vx)/2.d0 wy = (-AIMAG(b) - vy)/2.d0 ux = (-REAL(b) + vx)/2.d0 uy = (-AIMAG(b) + vy)/2.d0 r = REAL(a)*REAL(a) + AIMAG(a)*AIMAG(a) sx = (REAL(a)*wx + AIMAG(a)*wy)/r sy = (REAL(a)*wy - AIMAG(a)*wx)/r s1 = CMPLX(sx,sy) sx = (REAL(a)*ux + AIMAG(a)*uy)/r sy = (REAL(a)*uy - AIMAG(a)*ux)/r s2 = CMPLX(sx,sy) return End !solve complex equation a*z^3+b*z^2+c*z+d=0 and return !3 solutions z1,z2,z3. Subroutine Croot3(a,b,c,d, z1,z2,z3) Complex a,b,c,d,z1,z2,z3 !---------------------------------------------------- ! 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)) !------------------------------------------------------------------------------------ Complex p,q,det,sdet,u,u1,u2,v,v1,v2,w,zz Complex z(9) ! calculate p p = (c/a) - (b*b/(3*a*a)) ! calculate q q = (2*b**3/(27*a**3)) + (d/a) - (b*c/(3*a*a)) ! calculate det = q²/4 + p^3/27 det = (q*q/4.0) + (p*p*p/27.0) sdet = SQRT(det) !now sdet contains sqrt(q²/4 + p^3/27)) v = -q/2.0 + sdet call CUBRT(v,u,u1,u2) !3 cubic roots w = -q/2.0 - sdet call CUBRT(w,v,v1,v2) !3 cubic roots 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 index=1 do i=1, 9 ! call F1(p,q,z(i),zz) call F1(a,b,c,d,z(i),zz) if (CABS(zz) < 0.001) then !detect if solution is correct if (index==1) then z1=z(i); index=index+1 else if (index==2) then z2=z(i); index=index+1 else z3=z(i); end if end if end do return End ! solve complex equation a*z^4+b*z^3+c*z^2+d*z+e=0 and return ! 4 solutions z1,z2,z3,z4). Subroutine Croot4(a,b,c,d,e, z1,z2,z3,z4) Complex a,b,c,d,e,z1,z2,z3,z4 ! Labels: 100,200,300,400 Complex a1,aa,b1,bb,c1,cc, q,r,s,t,u,v,w Complex ck,l,m Complex z(3), zz(12) t=a; !save original a in t a=b/t; b=c/t; c=d/t; d=e/t ! Now equation is: z^4+a*z^3+b*z^2+c*z+d = 0 ! calculate q, r, s q = b - (3.0 * a * a / 8.0) r = c - (a * b / 2) + (a * a * a / 8.0) s = d - (a * c / 4) + (a * a * b / 16.0) - (3 * a * a * a * a / 256.0) ! Défine coefficients of cubic equation z^3+aa*z^2+bb*z+cc=0 (1) aa = q / 2.0 bb = (q * q - 4 * s) / 16.0 cc = -(r * r / 64.0) ! Calculate complex roots of equation (1) a1 = CMPLX(1.0, 0.0) IF (CABS(r) > 1E-15.OR.CABS(bb) > 1E-15) GOTO 100 ! Particular case when equation (1) is of 2nd order b1=aa; c1=bb call Equa2c(a1,b1,c1,z(1),z(2)) z(3) = CMPLX(0.0, 0.0) GOTO 200 100 call Croot3(a1,aa,bb,cc,z(1),z(2),z(3)) 200 do i=1, 3 ck = SQRT(z(i)) ! Calculate l and m if k=0 IF (CABS(ck) == 0.0) THEN r = SQRT(q * q - 4 * s) GOTO 300 end if q = q + (4 * z(i)) r = r / (2 * ck) 300 l = (q - r) / 2.0 m = (q + r) / 2.0 ! Solving two equations of degree 2 b1 = 2.0 * ck ! 1st equation call Equa2c(a1,b1,l,u,v) b1 =-b1 ! 2nd equation call Equa2c(a1,b1,m,w,t) ! Transferring solutions in zz(i) Select Case(i) Case(1) zz(1)=u zz(2)=v zz(3)=w zz(4)=t Case(2) zz(5)=u zz(6)=v zz(7)=w zz(8)=t Case(3) zz(9)=u zz(10)=v zz(11)=w zz(12)=t End Select !Note: only 4 solutions are correct end do !i loop do i=1, 12 !shift zz by -a/4 zz(i) = zz(i)-(a/4) end do ! select 4 correct solutions i=1 400 call F(a,b,c,d, zz(i),t) if (CABS(t) < 0.001) then Select Case(i) Case(1) z1=zz(1) z2=zz(2) z3=zz(3) z4=zz(4) Case(5) z1=zz(5) z2=zz(6) z3=zz(7) z4=zz(8) Case(9) z1=zz(9) z2=zz(10) z3=zz(11) z4=zz(12) End Select end if i = i + 4; if (i<=9) goto 400 ! now complex roots are in z1,z2,z3,z4. return End; ! end of file tcroot4.f90