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

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

[Fortran] syntaxhighlighter_viewsource syntaxhighlighter_copycode

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 00:07:48

问题有几个：
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] syntaxhighlighter_viewsource syntaxhighlighter_copycode

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

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

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

pasuka · 发表于 2014-3-20 16:44:45

正确的做法应该是转换为计算伴随矩阵的特征值
贴一个matlab代码：
function [x, w] = GaussLegendre_2(n)
% As you can see, the code consists of 2 blocks:
% 1: construct a symmetrical companion matrix
% 2: determine the (real) eigenvalues (i.e. the roots of the polynomial).
% It can produce the correct abscissas and weights, for any value n>=2.
%
% input: n ---- number of intergrate points
%
% output: x --- GaussLegendre intergrate point
% w --- GaussLegendre intergrate weight
%
i = 1:n-1;
a = i./sqrt(4*i.^2-1);
CM = diag(a,1) + diag(a,-1);
[V L] = eig(CM);
[x ind] = sort(diag(L));
V = V(:,ind)';
w = 2 * V(:,1).^2;
return
end

fcode · 发表于 2014-3-20 17:23:12

王培杰发表于 2014-3-20 15:40
勒让的多项式，N为奇数时，x=0必为零点，而module second的循环中m必定会过零点，只需把m初值设为-1.0000 ...

过零点的问题还是比较复杂的。
你得考虑浮点数的误差。如果恰好分段的节点附近存在零点就容易有问题。

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

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

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

[Fortran] syntaxhighlighter_viewsource syntaxhighlighter_copycode

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:23:58

我这里输出没问题，精度已经足够了。
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 就满足吧。

王培杰 · 发表于 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 21:28:16

把递推公式改为

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

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

新人勋章

爱心勋章

水王勋章

元老勋章

热心勋章

美女勋章

星光勋章

贡献勋章

管理勋章

帅哥勋章

规矩勋章