求勒让德多项式的零点，调试失败，求解哪里错了

王培杰 · 发表于 2014-3-21 22:42:03

chuxf 发表于 2014-3-21 21:28
把递推公式改为

fun(i)=((2*i-1)*x*fun(i-1)-(i-1)*fun(i-2))/i

终于成功了！谢谢@chuxf，帮我解决了很多问题。第二次编程终于完成了高斯勒让德积分法。一下是最后的完整代码，欢迎大家批评指正。

原理：高斯积分法 http://zh.wikipedia.org/wiki/高斯求积

实例：这里求sinx在（0，π/2）内的积分，结果应为1

程序如下：

[Fortran] 纯文本查看 复制代码

!高斯勒让德积分法
module first
    implicit none
    real,parameter::zero=1E-15
    integer,parameter::n=7
    contains
        real*8 function bisect(a,b)             !二分法精确求取勒让德多项式零点
        implicit none
        real*8::a,b,c,fa,fb,fc
        bisect=0d0
        do 
            c=(a+b)/2d0
            fa=func(a)
            fb=func(b)
            fc=func(c)
            if(fa*fc<0)then
                b=c
            else
                a=c
            end if
            if((b-a)<zero)exit
        end do 
        bisect=c
        end function

        real*8 function func(x)                 !求勒让德多项式的值
        implicit none
        real*8::x
        integer::i
        real*8::fun(n)
        fun(1)=x
        fun(2)=1.5*x*x-0.5
        do i=3,n
            fun(i)=((2*i-1)*x*fun(i-1)-(i-1)*fun(i-2))/i
        enddo
        func=fun(n)
        return
        end function
        
        real*8 function func1(x)                !求勒让德多项式的第n-1项值（不知道能不能并入上一个函数）
        implicit none
        real*8::x
        integer::i
        real*8::fun(n)
        fun(1)=x
        fun(2)=1.5*x*x-0.5
        do i=3,n
            fun(i)=((2*i-1)*x*fun(i-1)-(i-1)*fun(i-2))/i
        enddo
        func1=fun(n-1)
        return
        end function
        
        real*8 function fx(x)                   !被积函数f=sinx，（a，b）是积分区间，这里是（0，PI/2）
        implicit none
        real*8::x,y,a,b
        a=0
        b=1.57079632
        y=(a+b)/2+((b-a)/2)*x
        fx=sin(y)
        end function
end module  first

module second
use first
    implicit none
    contains
    subroutine fn0(fn)                          !对分法求勒让德多项式零点
        implicit none
        integer::i,j
        real*8 :: fn(:)
        real*8,allocatable :: k(:)
        real*8::p,q,m
        m=-1.0001_8
        j=1
        allocate(k(size(fn)))
        k=0.0_8
        do i=1,1999
            p=func(m)
            q=func(m+0.001_8)
            if(p*q<zero)then
                write(*,*)'j=',j,'m=',m 
                k(j)=m
                j=j+1
            endif
            m=m+0.001_8
        end do
        do i=1,j
            fn(i)=bisect(k(i),k(i)+0.001_8)     !调用二分法精确求解
        end do
    end subroutine
end module second
    
    

program GSLD                                    !高斯勒让德积分法
use first
use second
implicit none
integer::i
real*8::answer
real*8::pnn(n),ak(n),fn(n+1)
answer=0.0
call fn0(fn)
do i=1,n
    pnn(i)=(n*(func1(fn(i))-fn(i)*func(fn(i))))/(1-fn(i)**2d0)!求Pn的倒数
    ak(i)=2.0/n/func1(fn(i))/pnn(i)             !求Ak
    answer=answer+ak(i)*fx(fn(i))               !累加求和
    write(*,*)i,fn(i),ak(i)
end do
write(*,*)'answer=',answer*1.57079632/2
pause
endprogram GSLD

输出结果如图