问题有几个:
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] 纯文本查看 复制代码 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
|