论坛 › Fortran算法研究 › 求勒让德多项式的零点，调试失败，求解哪里错了

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

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

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 modulefirst

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
         58.812395257962180E-016
         60.396000000000001
         70.748000000000001
         80.960000000000002

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 modulefirst

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 = 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);
= eig(CM);
= 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.勒 ...

我把二分法也改好了，但输出的结果好像错了，请问是哪里出了问题。
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 modulefirst

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
         43.465581988049138E-016 -9.275231093379325E-016
         50.396520880908022      -1.966680786478849E-015
         60.747724940698209       1.332267629550188E-015
         70.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

查看完整版本: 求勒让德多项式的零点，调试失败，求解哪里错了