Studio di funzioni ellittiche
! ! Mio modo di calcolare le funzioni ellittiche di Jacobi. ! qui lavoro solo con dati reali. I valori con arg. ! complessi sono deducibili da quelli con arg. reali. ! ! I risultati vengono confrontati con la libreria IMSL che, ! mi sembra, AHIME', commettere degli errori quando il ! parametro m e' negativo. ! ! Il caso di m negativo e' piuttosto insolito e questo ! forse giustifica quello che ritengo un baco dell'IMSL. ! Oppure bisogna reputare che le formule dell' Abramowitz ! Stegun siano errate, il che sarebbe un fatto grave ! che incrinerebbe la credibilita' dell glorioso ! "Handbook of mathematical functions". ! ! Un indizio che sia l'IMSL a sbagliare e' il fatto che ! per m prossimo a zero, l'approssimazione in termini ! di funzioni circolari, per m negativo e' in disaccordo ! con quanto si ottiene con l'IMSL mentre e' in buon ! accordo con quanto si ottiene dalle formule ! dell' Abramowitz. ! ! Per u ed m piccoli infatti, approssimativamente: ! ! sn(u,m) = sin(u)-m*(u-sin(u)*cos(u))*cos(u)/4 ! cn(u,m) = cos(u)+m*(u-sin(u)*cos(u))*sin(u)/4 ! dn(u,m) = 1 - m*sin(u)*sin(u)/2 ! ! E le formule dell'Abramowitz concordano bene sia ! per m piccolo negativo che positivo.... ! ! Comunque ho realizzato anche una subroutine che ! calcola le funzioni ellittiche usando formule integrali ! canoniche ed anche essa conferma che l'errore e' nella ! libreria IMSL. ! ! Prendo l'occasione per evidenziare la formula che ho ! usato per integrare. Divido l'intervallino in 12 ! parti e calcolo la funzione integranda solo nei ! punti dispari. Poi sommo con pesi opportuni ossia ! uso la formula ! ! I = (247*f(1)+139*f(3)+254*f(5)+254*f(7)+ ! 139*f(9)+247*f(11))*delta/1280 ! ! Se la funzione e' regolare questa formula fornisce ! ottimi risultati numerici ed inoltre ha il pregio ! di non richiedere il calcolo della funzione ! integranda negli estremi dell'intervallo e pertanto ! possono essere integrate funzioni che vanno ! anche all'infinito negli estremi, purche' le ! funzioni siano veramente integrabili... ! module xscd implicit none contains subroutine agm(n,a,b,c,p) integer,intent(out)::n real(8),intent(in out),dimension(0:)::a,b,c real(8),optional,intent(in)::p real(8)::prec integer::j,ns ns=size(a,1) if(present(p)) then prec=abs(p) else prec=1.0d-14 end if passo: do j=0,ns-2 n=j+1 a(n)=a(j)+b(j) b(n)=sqrt(4.d0*a(j)*b(j)) c(n)=a(j)-b(j) if(prec>abs(c(n))) then exit passo end if end do passo return end subroutine agm ! recursive subroutine scd(snr,cnr,dnr,unr,pm) real(8),intent(in)::unr,pm real(8),intent(out)::snr,cnr,dnr real(8),dimension(0:100)::a,b,c real(8):: f,fh,se,cs,upm,pmm,sqm integer:: j,n if(1.d-9>abs(pm-1.d0)) then f=1.0d0-pm se=sinh(unr) cs=cosh(unr) snr= se/cs +0.25d0*f*(se*cs-unr)/(cs*cs) cnr=1.d0/cs-0.25d0*f*(se*cs-unr)*se/(cs*cs) dnr=1.d0/cs+0.25d0*f*(se*cs+unr)*se/(cs*cs) return else if(1.d-10>abs(pm)) then se=sin(unr) cs=cos(unr) snr=se-0.25d0*pm*(unr-se*cs)*cs cnr=cs+0.25d0*pm*(unr-se*cs)*se dnr=1.d0-0.5d0*pm*se*se return else if(pm>1.0d0) then upm=1.0/pm sqm=sqrt(pm) call scd (snr,cnr,dnr,unr*sqm,upm) snr=snr/sqm f=dnr dnr=cnr cnr=f return else if (0.d0>pm) then pmm=-pm upm=1.0d0/(1.0d0+pmm) sqm=sqrt(upm) call scd (snr,cnr,dnr,unr/sqm,pmm*upm) dnr=1.0d0/dnr cnr=cnr*dnr snr=snr*dnr*sqm return end if a(0)=1.0d0 b(0)=sqrt(a(0)-pm) c(0)=sqrt(pm) call agm(n,a,b,c) f=a(n)*unr if(n > 10) then write(*,*) "n=",n end if do j=n,1,-1 fh=f f=0.5d0*(asin(c(j)*sin(f)/a(j))+f) end do snr=sin(f) cnr=cos(f) dnr=cnr/cos(fh-f) return end subroutine scd end module xscd ! module integro implicit none ! ! Applica il modo INTEGRALE per calcolare le ! funzioni ellittiche. ! contains ! subroutine trova_u(s,c,d,u,fi,pm,nstep) ! ! s=sn(u,pm) ; c=cn(u,pm) ; d=dn(u,pm) ! ! ma u, a sua volta non e' indipendente ma ! dipende da fi. ! real(8),intent(out)::s,c,d,u integer,intent(in):: nstep real(8),intent(in):: fi,pm real(8):: delta,delta6,delta12,t integer::j if(0.d0>=abs(fi)) then u=0.d0 return end if delta=fi/nstep u=0.d0 delta12=delta/12.d0 delta6=2.0d0*delta12 ! ! Passi di integrazione... ! do j=0,nstep-1 t=j*delta +delta12 u= u+ 247.0d0*delta/sqrt(1.d0-pm*sin(t)*sin(t)) t=t+delta6 u= u+ 139.0d0*delta/sqrt(1.d0-pm*sin(t)*sin(t)) t=t+delta6 u= u+ 254.0d0*delta/sqrt(1.d0-pm*sin(t)*sin(t)) t=t+delta6 u= u+ 254.0d0*delta/sqrt(1.d0-pm*sin(t)*sin(t)) t=t+delta6 u= u+ 139.0d0*delta/sqrt(1.d0-pm*sin(t)*sin(t)) t=t+delta6 u= u+ 247.0d0*delta/sqrt(1.d0-pm*sin(t)*sin(t)) end do u=u/1280.0d0 s=sin(fi) c=cos(fi) d=sqrt(1-pm*s*s) return end subroutine trova_u ! end module integro ! program testscd use msimsl use xscd use integro implicit none real(8):: snr,cnr,dnr,unr,pm,fi real(8):: s,c,d,vmax=1.0d-2 real(8)::imsl_sn,imsl_cn,imsl_dn integer::j,k write(unit=*,fmt=*) & "Verifica delle formule dell'Abramowitz-Stegun sulle" write(unit=*,fmt=*) & "funzioni ellittiche jacobiane." write(unit=*,fmt=*) & "Viene usata anche la libreria IMSL che pero' risulta" write(unit=*,fmt=*) & "errata !!!! quando il valore del parameto pm e' negativo" write(unit=*,fmt=*) " " studio:do write(unit=*,fmt=*) & "Assegna pm (parametro delle sn(u,pm) etc...: 0==fine) " read(unit=*,fmt=*) pm if(0.d0>= abs(pm) ) then exit studio end if write(unit=*,fmt=*) "Lavora con",pm do j=1,20 fi=j*0.1d0 if(1.0d0 -pm > vmax) then call trova_u(s,c,d,unr,fi,pm,50) else unr=fi d=0.0 end if call scd (snr,cnr,dnr,unr,pm) imsl_sn=dejsn(unr,pm) imsl_cn=dejcn(unr,pm) imsl_dn=dejdn(unr,pm) write(unit=*,fmt="(1x,g12.5,a,g12.5,a)") & unr," == u ; ",dnr-dejdn(unr,pm)," == diff. da IMSL" if(1.0d0 - pm > vmax) then write(unit=*,fmt="(1x,g12.5,a,g12.5,a,g12.5,a)") & snr," == sn(u); ",s," == test; ",snr-s," == sn(u)-test " end if end do end do studio write(unit=*,fmt=*) "Le funzioni dejsn(x,m), dejcn(x,m), dejdn(x,m)" write(unit=*,fmt=*) "della libreria IMSL acclusa al Digital Visual" write(unit=*,fmt=*) "Fortran. versione 5, per m negativo.... sono" write(unit=*,fmt=*) "dunque bacate !!!!! :-( " write(unit=*,fmt=*) " " write(unit=*,fmt=*) "Amen. Invia...." read(unit=*,fmt=*) stop end program testscd ! !