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

楚香饭

问题有几个：
1.可分配数组，必须要经过分配然后才能使用。并不是说可分配数组就是可随意伸缩的数组。
2.勒让德多项式的递推公式写错了。少写一个 x，而且顺序写错了，注意乘法和除法的优先顺序。
(n+1)P_{n+1}=(2n+1)xP_{n}-nP_{n-1}
3.fn0 函数用来寻找0点附近的区域，在0的位置会产生两个相同的“区间”，实际0点只有一个。这个要自己想办法来解决。
4.后面的二分法求零点，貌似也有点问题。太晚了，明天再调试。

以下程序我自己修改的。能得出 8 个区间（实际7个零点）。错误我都标在注释里了，供您参考

输出的结果是：
           1 -0.961000000000000
           2 -0.749000000000000
           3 -0.396999999999999
           4 -9.999999999991188E-004
           5  8.812395257962180E-016
           6  0.396000000000001
           7  0.748000000000001
           8  0.960000000000002

[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
        bisect = 0.0
        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)*x*fun(i-1)-n*fun(i-2) ) / (n+1)  !// 递推公式错误，少写一个 x ，注意 n+1 要最后除
        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(kind=8)::p,q,m
        m=-1.0_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,m !// 把变号区间输出
                k(j)=m
                j=j+1
            endif
            m=m+0.001_8
        end do
        do i=1,j-1
            fn(i)=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(:)
allocate(fn(2000)) !// fn 需要先分配
call fn0(fn)
!write(*,*)fn(1:100) !// 暂时不输出fn
pause
endprogram GSLD

楚香饭

我这里输出没问题，精度已经足够了。
j=           1 m= -0.959100000000000
j=           2 m= -0.748100000000000
j=           3 m= -0.397099999999999
j=           4 m= -9.999999999910775E-005
j=           5 m=  0.395900000000001
j=           6 m=  0.746900000000001
j=           7 m=  0.957900000000002
           1 -0.958533820646988      -2.791417890486108E-015
           2 -0.747724940698208      -3.616154994493367E-015
           3 -0.396520880908023      -1.395708945243054E-015
           4  3.465581988049138E-016 -9.275231093379325E-016
           5  0.396520880908022      -1.966680786478849E-015
           6  0.747724940698209       1.332267629550188E-015
           7  0.958533820646989       1.034093445793717E-014

红色是反算的值，非常接近 0
浮点数是有误差的，所以精确到 e-14,-15 就满足吧。

楚香饭

把递推公式改为

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

楚香饭

pasuka 发表于 2014-3-27 21:49
lz可以尝试把多项式阶次提升到50以上，若不采用计算特征值的办法，计算出来的积分点和权系数是不正确的 ...

从算法角度上来说，楼主的代码确实有局限性。pasuka 的算法更合理。

根据楼主的代码看，貌似只是学习。