Fortran Coder

标题: 求勒让德多项式的零点，调试失败，求解哪里错了 [打印本页]

---

作者: 王培杰    时间: 2014-3-19 22:48
标题: 求勒让德多项式的零点，调试失败，求解哪里错了

[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

QQ图片20140319224616.jpg (56.93 KB, 下载次数: 439)

下载附件

2014-3-19 22:47 上传

---

作者: 楚香饭    时间: 2014-3-20 00:07
问题有几个：
1.可分配数组，必须要经过分配然后才能使用。并不是说可分配数组就是可随意伸缩的数组。
2.勒让德多项式的递推公式写错了。少写一个 x，而且顺序写错了，注意乘法和除法的优先顺序。
[latex](n+1)P_{n+1}=(2n+1)xP_{n}-nP_{n-1}[/latex]
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

---

作者: 王培杰    时间: 2014-3-20 15:40

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

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

QQ图片20140320153925.jpg (23.98 KB, 下载次数: 407)

下载附件

2014-3-20 15:40 上传

---

作者: 王培杰    时间: 2014-3-20 15:44

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

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

---

作者: pasuka    时间: 2014-3-20 16:44
正确的做法应该是转换为计算伴随矩阵的特征值
贴一个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

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

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

---

作者: 王培杰    时间: 2014-3-21 20:11

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

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

结果和书上的出入很大

QQ图片20140321202849.jpg (97.7 KB, 下载次数: 425)

下载附件

2014-3-21 20:29 上传

---

作者: 楚香饭    时间: 2014-3-21 21:28
把递推公式改为

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

---

作者: 王培杰    时间: 2014-3-21 22:42

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

输出结果如图

QQ图片20140321222927.jpg (63 KB, 下载次数: 336)

下载附件

2014-3-21 22:41 上传

---

作者: aliouying    时间: 2014-3-27 21:31
明天我测试下我程序中的勒让德高斯积分是否精度足够

---

作者: fcode    时间: 2014-3-27 21:38
楼上赶着审核资料还来论坛，敬业精神可歌可泣

---

作者: pasuka    时间: 2014-3-27 21:49

王培杰 发表于 2014-3-21 22:42
终于成功了！谢谢@chuxf，帮我解决了很多问题。第二次编程终于完成了高斯勒让德积分法。一下是最后的完整 ...

lz可以尝试把多项式阶次提升到50以上，若不采用计算特征值的办法，计算出来的积分点和权系数是不正确的

---

作者: 楚香饭    时间: 2014-3-27 21:58

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

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

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

---

欢迎光临 Fortran Coder (http://bbs.fcode.cn/)

Powered by Discuz! X3.2