Fortran傅里叶级数逼近

2020-08-30 本文已影响0人 zoziha

1 傅里叶级数逼近

1.1 Fortran源码

program main

    dimension f(81), a(41), b(41)
    double precision f,a,b,c,h

    write(*,*)
    write(*,10)
    10  format(3x,'n',12x,'a(n)',16x,'b(n)')
    h=8*atan(1.)/81.0
    do i=1,81
      c=h*(i-0.5)
      f(i)=c*c
    end do
    n1=41
    n2=81
    call kfour(f,n1,n2,a,b)
    do i=1,n1
      k=i-1
      write(*,40) k,a(i),b(i)
    end do
    write(*,*)
    40  format(1x,i3,5x,d15.6,4x,d15.6)
    call firm(f,a,b)

end program
subroutine kfour(f,n1,n2,a,b)

    dimension f(n2),a(n1),b(n1)
    double precision f,a,b,c,s,t,c1,s1,u1,u2,u0
    real*8,parameter :: PI=4*atan(1.)

    t=2*PI/n2
    c=cos(t)
    s=sin(t)
    t=2.0/n2
    c1=1.0
    s1=0.0
        do i=1,n1
            u1=0.0
            u2=0.0
            do j=n2,2,-1
                u0=f(j)+2.0*c1*u1-u2
                u2=u1
                u1=u0
            end do
            a(i)=t*(f(1)+u1*c1-u2)
            b(i)=t*u1*s1
            u0=c*c1-s*s1
            s1=c*s1+s*c1
            c1=u0
        end do
    return

end subroutine
subroutine firm(f,a,b)

    real*8 :: f(81),a(41),b(41),fo(81)

    h=8*atan(1.)/81.0
    fo=0
    do i=1,81
        c=h*(i-0.5)
        do j=2,41
            fo(i)=fo(i)+a(j)*cos((j-1)*c)+b(j)*sin((j-1)*c)
        end do
        fo(i)=fo(i)+0.5*a(1)
    end do

    do i=1,81
        write(*,100) h*(i-0.5),fo(i),f(i)
100     format(3(ES12.3E3,2X))
    end do

end subroutine