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

王培杰

[Fortran] 纯文本查看 复制代码

查看代码

打印?

01	module first

02	implicit none

03	real,parameter::zero=1E-8

04	integer,parameter::n=7

05

contains

06	real function bisect(a,b)

07	implicit none

08	real*8::a,b,c,fa,fb,fc

09	c=(a+b)/2.0

10	fc=func(c)

11	do while(abs(fc)>zero)

12	fa=func(a)

13	fb=func(b)

14	if(fa*fb<0)then

15

b=c

16	c=(a+b)/2.0

17

else

18

a=c

19	c=(a+b)/2.0

20

end if

21	fc=func(c)

22

return

23

end do

24

bisect=c

25

return

26	end function

27

28	real function func(x)

29	implicit none

30

real*8::x

31	integer::i

32	real*8::fun(n)

33

fun(1)=x

34	fun(2)=1.5*x*x-0.5

35

do i=3,n

36	fun(i)=(2*n+1)/(n+1)*fun(i-1)-n/(n+1)*fun(i-2)

37

enddo

38	func=fun(n)

39

return

40	end function

41	end module  first

42

43	module second

44

use first

45	implicit none

46

contains

47	subroutine fn0(fn)

48	implicit none

49	integer::i,j

50	real*8,allocatable::fn(:),k(:)

51	real*8::p,q,m

52

m=-1

53

j=1

54	do i=1,1999

55

p=func(m)

56	q=func(m+0.001)

57	if(p*q<zero)then

58

k(j)=m

59

j=j+1

60

endif

61

m=m+0.001

62

end do

63	do i=1,j-1

64	fn(1)=bisect(k(j),k(j)+0.001)

65

end do

66	end subroutine

67	end module second

68

69

70	program GSLD

71

use first

72	use second

73	implicit none

74	real*8,allocatable::fn(:)

75	call fn0(fn)

76	write(*,*)fn

77

pause

78	endprogram GSLD

楚香饭

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

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

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

pasuka

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

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

fcode

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

aliouying

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

王培杰

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] 纯文本查看 复制代码

查看代码

打印?

001

!高斯勒让德积分法

002	module first

003	implicit none

004	real,parameter::zero=1E-15

005	integer,parameter::n=7

006

contains

007	real*8 function bisect(a,b)             !二分法精确求取勒让德多项式零点

008	implicit none

009	real*8::a,b,c,fa,fb,fc

010	bisect=0d0

011

do

012	c=(a+b)/2d0

013	fa=func(a)

014	fb=func(b)

015	fc=func(c)

016	if(fa*fc<0)then

017

b=c

018

else

019

a=c

020

end if

021	if((b-a)<zero)exit

022

end do

023

bisect=c

024	end function

025

026	real*8 function func(x)                 !求勒让德多项式的值

027	implicit none

028

real*8::x

029	integer::i

030	real*8::fun(n)

031

fun(1)=x

032	fun(2)=1.5*x*x-0.5

033

do i=3,n

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

035

enddo

036	func=fun(n)

037

return

038	end function

039

040	real*8 function func1(x)                !求勒让德多项式的第n-1项值（不知道能不能并入上一个函数）

041	implicit none

042

real*8::x

043	integer::i

044	real*8::fun(n)

045

fun(1)=x

046	fun(2)=1.5*x*x-0.5

047

do i=3,n

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

049

enddo

050	func1=fun(n-1)

051

return

052	end function

053

054	real*8 function fx(x)                   !被积函数f=sinx，（a，b）是积分区间，这里是（0，PI/2）

055	implicit none

056	real*8::x,y,a,b

057

a=0

058	b=1.57079632

059	y=(a+b)/2+((b-a)/2)*x

060

fx=sin(y)

061	end function

062	end module  first

063

064	module second

065

use first

066	implicit none

067

contains

068	subroutine fn0(fn)                          !对分法求勒让德多项式零点

069	implicit none

070	integer::i,j

071	real*8 :: fn(:)

072	real*8,allocatable :: k(:)

073	real*8::p,q,m

074	m=-1.0001_8

075

j=1

076	allocate(k(size(fn)))

077

k=0.0_8

078	do i=1,1999

079

p=func(m)

080	q=func(m+0.001_8)

081	if(p*q<zero)then

082	write(*,*)'j=',j,'m=',m

083

k(j)=m

084

j=j+1

085

endif

086	m=m+0.001_8

087

end do

088

do i=1,j

089	fn(i)=bisect(k(i),k(i)+0.001_8)     !调用二分法精确求解

090

end do

091	end subroutine

092	end module second

093

094

095

096	program GSLD                                    !高斯勒让德积分法

097

use first

098	use second

099	implicit none

100	integer::i

101	real*8::answer

102	real*8::pnn(n),ak(n),fn(n+1)

103	answer=0.0

104	call fn0(fn)

105

do i=1,n

106	pnn(i)=(n*(func1(fn(i))-fn(i)*func(fn(i))))/(1-fn(i)**2d0)!求Pn的倒数

107	ak(i)=2.0/n/func1(fn(i))/pnn(i)             !求Ak

108	answer=answer+ak(i)*fx(fn(i))               !累加求和

109	write(*,*)i,fn(i),ak(i)

110

end do

111	write(*,*)'answer=',answer*1.57079632/2

112

pause

113	endprogram GSLD

输出结果如图

楚香饭

把递推公式改为

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

王培杰

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

结果和书上的出入很大

楚香饭

我这里输出没问题，精度已经足够了。
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] 纯文本查看 复制代码

查看代码

打印?

01	module first

02	implicit none

03	real,parameter::zero=1E-15

04	integer,parameter::n=7

05

contains

06	real*8 function bisect(a,b)

07	implicit none

08	real*8::a,b,c,fa,fb,fc

09	bisect=0d0

10

do

11	c=(a+b)/2d0

12	fa=func(a)

13	fb=func(b)

14	fc=func(c)

15	if(fa*fc<0)then

16

b=c

17

else

18

a=c

19

end if

20	if((b-a)<zero)exit

21

end do

22

bisect=c

23	end function

24

25	real*8 function func(x)

26	implicit none

27

real*8::x

28	integer::i

29	real*8::fun(n)

30

fun(1)=x

31	fun(2)=1.5*x*x-0.5

32

do i=3,n

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

34

enddo

35	func=fun(n)

36

return

37	end function

38	end module  first

39

40	module second

41

use first

42	implicit none

43

contains

44	subroutine fn0(fn)

45	implicit none

46	integer::i,j

47	real*8 :: fn(:)

48	real*8,allocatable :: k(:)

49	real*8::p,q,m

50	m=-1.0001_8

51

j=1

52	allocate(k(size(fn)))

53

k=0.0_8

54	do i=1,1999

55

p=func(m)

56	q=func(m+0.001_8)

57	if(p*q<zero)then

58	write(*,*)'j=',j,'m=',m

59

k(j)=m

60

j=j+1

61

endif

62	m=m+0.001_8

63

end do

64

do i=1,j

65	fn(i)=bisect(k(i),k(i)+0.001_8)

66

end do

67	end subroutine

68	end module second

69

70

71	program GSLD

72

use first

73	use second

74	implicit none

75	integer::i

76	real*8,allocatable::fn(:)

77	allocate(fn(200))

78	call fn0(fn)

79

do i=1,n

80	write(*,*)i,fn(i)

81

enddo

82

pause

83	endprogram GSLD