把递推公式改为
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
程序如下：
!高斯勒让德积分法
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 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::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
输出结果如图
明天我测试下我程序中的勒让德高斯积分是否精度足够 楼上赶着审核资料还来论坛，敬业精神可歌可泣 王培杰 发表于 2014-3-21 22:42
终于成功了！谢谢@chuxf，帮我解决了很多问题。第二次编程终于完成了高斯勒让德积分法。一下是最后的完整 ...
lz可以尝试把多项式阶次提升到50以上，若不采用计算特征值的办法，计算出来的积分点和权系数是不正确的 pasuka 发表于 2014-3-27 21:49
lz可以尝试把多项式阶次提升到50以上，若不采用计算特征值的办法，计算出来的积分点和权系数是不正确的 ...
从算法角度上来说，楼主的代码确实有局限性。pasuka 的算法更合理。
根据楼主的代码看，貌似只是学习。
页:
1
[2]