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

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

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

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

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

楼上赶着审核资料还来论坛，敬业精神可歌可泣

明天我测试下我程序中的勒让德高斯积分是否精度足够

把递推公式改为

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

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

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

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

[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