module GSLD !高斯—勒让德积分高斯点及权重的求解模块
implicit none
integer,parameter::n=5 !设置求解高斯点的个数
integer,parameter::nn=16 !设置kind的数值
contains
real(kind=nn) function f(x) !定义四阶精度的函数f（x）
implicit none
integer::i
real(kind=nn)::a(n),x !a(n)代表n阶勒让德多项式
a(1)=x !1阶勒让德多项式
a(2)=1.5_16*(x**2)-0.5_16 !2阶勒让德多项式
do i=3,n
a(i)=(2*i-1)*x*a(i-1)/i-(i-1)*a(i-2)/i
!利用递推关系产生n阶勒让德多项式
enddo
f=a(n) !生成的n阶勒让德多项式
end function
real(kind=nn) function f1(x)
implicit none
integer::i
real(kind=nn)::a(n),x
a(1)=x
a(2)=1.5_16*x**2-0.5_16
do i=3,n-1
a(i)=(2*i-1)*x*a(i-1)/i-(i-1)*a(i-2)/i
enddo
f1=a(n-1) !生成的（n-1）阶勒让德多项式
endfunction
real(kind=nn) function g(x)
implicit none
integer::i
real(kind=nn)::a(n),x
a(1)=x
a(2)=1.5_16*x**2-0.5_16
do i=3,n
a(i)=(2*i-1)*x*a(i-1)/i-(i-1)*a(i-2)/i
enddo
g=n*a(n-1)/(1-x**2)-n*x*a(n)/(1-x**2)
!生成n阶勒让德多项式的导数表达式
endfunction
real(kind=nn) function bis(a,b) !二分法求解函数的解
implicit none
real(kind=nn)::a,b,c
!a,b是传递进来的划分好的有一个解存在的区间
do while(.true.)
c=(a+b)/2.0_16
if(f(c)*f(a)<0)then
b=c
else
a=c
endif
if((b-a)<1E-34)exit !这个精度设置的。。（有点狠了）
enddo
bis=c !bis即是利用二分法求得的解
end function
subroutine fn0(fn,ak)
implicit none
real(kind=nn)::m,fn(n),ak(n)
!定义数组,大小n由module开始声明。
integer::i,j
open(10,file='高斯点+权重.dat')
j=0 !赋值控制循环变量的初值
m=-1.000001 !设置计算域[-1，1] 的下限，即代替-1
do i=1,200000 !这个循环次数应该是由步长0.00001决 定,计算方法：200000=2/0.000001
if(f(m)*f(m+0.00001)<0)then !从下限处开始往上逐步累加，
!由步长0.00001说明最多求解10^5个解
j=j+1 !记录这是第几个解
fn(j)=bis(m,m+0.00001)
!调用二分法求解程序在分好的一小段上求解，
!将解存储在fn（j）
ak(j)=2.0_16/(n*f1(fn(j))*g(fn(j))) !高斯点的权重
write(*,*) '高斯点序号',j
write(*,*) '高斯点' ,fn(j)
write(*,*) '高斯点权重',ak(j)
write(10,*) '高斯点',fn(j)
write(10,*) '权重' ,ak(j)
endif
m=m+0.00001 !执行完一次判断m向前推进一步
enddo
close (10)
endsubroutine
real(kind=nn) function y(x) !被积函数
implicit none
real(kind=nn)::x
y=x**2/3.0_16 !每次计算在这里修改被积函数即可
end function
end module
program www_fcode_cn
use GSLD
implicit none
real(kind=nn)::fn(n),ak(n),x,a,b,answer
integer::i
open(11,file='GS_LD_result.dat')
call fn0(fn,ak) !调用求高斯零点和权重的子函数
a=0.0_16 !积分上限
b=1.0_16 !积分下限
answer=0.0_16 !求积分结果赋初始值
do i=1,n
answer=answer+ak(i)*y((a+b)/2.0_16+(b-a)/2.0_16*fn(i))
!部分高斯勒让德求积分公式
enddo
answer=answer*(b-a)/2.0_16
write(11,*) '高斯_勒让德求积分结果',answer
write(*,*) '高斯_勒让德求积分结果',answer
close(11)
stop
end program www_fcode_cn
module GSLD
implicit none
integer,parameter::n=5
contains
real*16 function f(x)
integer::i
real*16::a(n),x
a(1)=x
a(2)=1.5*(x**2)-0.5
do i=3,n
a(i)=(2*i-1)*x*a(i-1)/i-(i-1)*a(i-2)/i
enddo
f=a(n)
endfunction
real*16 function f1(x)
integer::i
real*16::a(n),x
a(1)=x
a(2)=1.5*x**2-0.5
do i=3,n
a(i)=(2*i-1)*x*a(i-1)/i-(i-1)*a(i-2)/i
enddo
f1=a(n-1)
endfunction
real*16 function g(x)
integer::i
real*16::a(n),x
a(1)=x
a(2)=1.5d0*x**2-0.5d0
do i=3,n
a(i)=(2*i-1)*x*a(i-1)/i-(i-1)*a(i-2)/i
enddo
g=n*a(n-1)/(1-x**2)-n*x*a(n)/(1-x**2)
endfunction
real*16 function bis(a,b)
real*16::a,b,c
do
c=(a+b)/2d0
if(f(c)*f(a)<0)then
b=c
else
a=c
endif
if((b-a)<1E-34)exit
enddo
bis=c
endfunction
subroutine fn0(fn,ak)
real*16::m,fn(n),ak(n)
integer::i,j
j=0
m=-1.00000001
do i=1,200000
if(f(m)*f(m+0.00001d0)<0)then
j=j+1
fn(j)=bis(m,m+0.00001)
ak(j)=2/(n*f1(fn(j))*g(fn(j)))
write(*,*)j,fn(j),ak(j)
endif
m=m+0.00001
enddo
endsubroutine
real*16 function y(x,x1,x2) !被积函数
real*16::x,x1,x2
y=2*x+2*x1+2*x2
endfunction
endmodule
program www_fcode_cn
use GSLD
implicit none
real*16::fn(n),ak(n),x,x1,x2,a,b,a1,b1,a2,b2,answer
integer::i,j,k
call fn0(fn,ak)
a=0
b=1 !积分上下限
a1=0
b1=1
a2=0
b2=1
answer=0
do k=1,n
do j=1,n
do i=1,n answer=answer+ak(i)*ak(j)*ak(k)*y((a+b)/2+(b-a)/2*fn(i),(a1+b1)/2+(b1-a1)/2*fn(j),(a1+b1)/2+(b1-a1)/2*fn(k))
enddo
enddo
enddo
answer=answer*(b-a)/2*(b1-a1)/2*(b2-a2)/2
write(*,*)answer
pause
end program www_fcode_cn
!x(n)高斯零点;w(n)权系数;n高斯点个数(可变)
program main
implicit none
integer,parameter:: n=1000 !n高斯点个数(可变)
integer,parameter::nn=8
real(kind=nn) ::x1,x2,x(n),w(n)
real(kind=nn) ::y,f,a
integer ::i
!------------------------------------------------
!x1,x2积分上下限;变换积分上下限为[-1,1 ]
x1=-10.0e1
x2=10.0e1
call gauleg(-1.d0,1.d0,x,w,n) !调用求解高斯点程序
y=0.0 !积分结果赋初始值
do i=1,n
y=y+0.5*(x2-x1)*w(i)*f(0.5*(x1+x2)+0.5*(x2-x1)*x(i)) !高斯求积公式，已经对被积函数做变量代换处理，
!将被积分区间换到【-1,1】
end do
write(*,*) '积分理论值',sqrt(3.1415926)*exp(-1.d0/4.d0) !测试计算理论值
write(*,*) '数值计算值',y !测试数值计算值
end
function f(x)
implicit none
integer,parameter::nn=8
real(kind=nn):: f,x
f=exp(-x**2)*cos(x)
end function
SUBROUTINE gauleg(x1,x2,x,w,n) !gauleg(-1.0,1.0,x,w,n)
!这x1，x2的传不传的没意思吧，勒让德求0点不就是在[-1,1]上吗？
implicit none
integer::n !计算的阶数
integer,parameter::nn=8
real(kind=nn)::x1,x2,x(n),w(n)!积分上下限传递进来
real(kind=nn),parameter::EPS=3.d-14!系数为何取3，真心诡异
integer::i,j,m
real(kind=8)::p1,p2,p3,pp,xl,xm,z,z1
m=(n+1)/2 !截取一半，这里取整，这个很厉害
xm=0.5d0*(x2+x1)
xl=0.5d0*(x2-x1) !将求解区间截取一半，在一半区间上求解高斯点
do 12 i=1,m !只计算一半，这个牛逼[0,1]区间求解
z=cos(3.141592654d0*(i-0.25d0)/(n+0.5d0)) !这是尼玛到底是啥？？
1 continue
p1=1.0
p2=0.0
do 11 j=1,n
p3=p2
p2=p1
p1=((2.0*j-1.0)*z*p2-(j-1.0)*p3)/j
11 continue
pp=n*(z*p1-p2)/(z*z-1.0) !代替n阶勒让德多项式的导数？？？
z1=z
z=z1-p1/pp
if(abs(z-z1).gt.EPS)goto 1
x(i)=xm-xl*z !代替高斯点
x(n+1-i)=xm+xl*z !代替另外一半的高斯点
w(i)=2.d0*xl/((1.d0-z*z)*pp*pp) !权重的表达式
w(n+1-i)=w(i) !对称的另一半的权重
12 continue
print*, x, w
return
END
QQ图片20140617120026.jpg (13.67 KB, 下载次数: 685)
pasuka 发表于 2014-6-16 18:08
lz文献查阅力度还不够，对于振荡函数的积分请用自适应高斯-克罗朗德积分一般来说21个点的高斯-克罗朗德积分 ...
岸边的鱼 发表于 2014-6-17 19:41
终于看懂你的lz是啥意思了，原来是楼主啊，自适应法不知道是不是就是自动变步长法？我其实是有使用的，不 ...
pasuka 发表于 2014-6-18 09:41
抱歉，单位网络没法上qq
至于gauss-kronrod的基本原理，lz可以自己网上搜索或者查看matlab的quadgk函数代 ...
瑶远梦想 发表于 2014-8-3 22:18
个人觉得用迭代改变步长时解决震荡函数积分的好方法！
| 欢迎光临 Fortran Coder (http://bbs.fcode.cn/) | Powered by Discuz! X3.2 |
