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

王培杰 · 发表于 2014-3-19 22:48:26

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

module first
    implicit none
    real,parameter::zero=1E-8
    integer,parameter::n=7
    contains
        real function bisect(a,b)
        implicit none
        real*8::a,b,c,fa,fb,fc
        c=(a+b)/2.0
        fc=func(c)
        do while(abs(fc)>zero)
            fa=func(a)
            fb=func(b)
            if(fa*fb<0)then
                b=c
                c=(a+b)/2.0
            else
                a=c
                c=(a+b)/2.0
            end if
            fc=func(c)
            return
        end do 
        bisect=c
        return
        end function

        real 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*n+1)/(n+1)*fun(i-1)-n/(n+1)*fun(i-2)
        enddo
        func=fun(n)
        return
        end function
end module  first

module second
use first
    implicit none
    contains
    subroutine fn0(fn)
        implicit none
        integer::i,j
        real*8,allocatable::fn(:),k(:)
        real*8::p,q,m
        m=-1
        j=1
        do i=1,1999
            p=func(m)
            q=func(m+0.001)
            if(p*q<zero)then
                k(j)=m
                j=j+1
            endif
            m=m+0.001
        end do
        do i=1,j-1
            fn(1)=bisect(k(j),k(j)+0.001)
        end do
    end subroutine
end module second
    
    
program GSLD
use first
use second
implicit none
real*8,allocatable::fn(:)
call fn0(fn)
write(*,*)fn
pause
endprogram GSLD

王培杰 · 发表于 2014-3-20 15:40:43

chuxf 发表于 2014-3-20 00:07
问题有几个：
1.可分配数组，必须要经过分配然后才能使用。并不是说可分配数组就是可随意伸缩的数组。
2.勒 ...

勒让的多项式，N为奇数时，x=0必为零点，而module second的循环中m必定会过零点，只需把m初值设为-1.00001（使其在循环中不过零点）即可。谢谢了！

王培杰 · 发表于 2014-3-20 15:44:18

chuxf 发表于 2014-3-20 00:07
问题有几个：
1.可分配数组，必须要经过分配然后才能使用。并不是说可分配数组就是可随意伸缩的数组。
2.勒 ...

还有一点请教一下，怎么能尽量提高这个程序结果的精度？

王培杰 · 发表于 2014-3-21 20:11:35

chuxf 发表于 2014-3-20 00:07
问题有几个：
1.可分配数组，必须要经过分配然后才能使用。并不是说可分配数组就是可随意伸缩的数组。
2.勒 ...

我把二分法也改好了，但输出的结果好像错了，请问是哪里出了问题。

[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*n-1)*x*fun(i-1)-(n-1)*fun(i-2))/n
        enddo
        func=fun(n)
        return
    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,allocatable::fn(:)
allocate(fn(200)) 
call fn0(fn)
do i=1,n
    write(*,*)i,fn(i) 
enddo
pause
endprogram GSLD

王培杰 · 发表于 2014-3-21 20:30:02

chuxf 发表于 2014-3-21 20:23
我这里输出没问题，精度已经足够了。
j= 1 m= -0.959100000000000
j= 2 m= -0.74810 ...

结果和书上的出入很大

王培杰 · 发表于 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

输出结果如图

[特殊函数] 求勒让德多项式的零点，调试失败，求解哪里错了