[Fortran] 纯文本查看 复制代码
Module gsld !高斯—勒让德积分高斯点及权重的求解模块
Implicit None
Integer, Parameter :: n = 5 !设置求解高斯点的个数
Integer, Parameter :: DP = Selected_Real_Kind( p=13 ) !设置kind的数值
Real(kind=DP) , parameter :: EPS = 1.0e-15_DP !精度设置
Contains
Real(Kind=DP) Function f(x) !定义函数f（x）
Implicit None
Integer :: i
Real (Kind=DP) :: a(n), x !a(n)代表n阶勒让德多项式
a(1) = x !1阶勒让德多项式
a(2) = 1.5_DP*(x**2) - 0.5_DP !2阶勒让德多项式
Do i = 3, n
a(i) = (2*i-1)*x*a(i-1)/i - (i-1)*a(i-2)/i
!利用递推关系产生n阶勒让德多项式
End Do
f = a(n) !生成的n阶勒让德多项式
End Function f
Real(Kind=DP) Function f1(x)
Implicit None
Integer :: i
Real (Kind=DP) :: a(n), x
a(1) = x
a(2) = 1.5_DP*x**2 - 0.5_DP
Do i = 3, n - 1
a(i) = (2*i-1)*x*a(i-1)/i - (i-1)*a(i-2)/i
End Do
f1 = a(n-1) !生成的（n-1）阶勒让德多项式
End Function f1
Real (Kind=DP) Function g(x)
Implicit None
Integer :: i
Real (Kind=DP) :: a(n), x
a(1) = x
a(2) = 1.5_DP*x**2 - 0.5_DP
Do i = 3, n
a(i) = (2*i-1)*x*a(i-1)/i - (i-1)*a(i-2)/i
End Do
g = n*a(n-1)/(1-x**2) - n*x*a(n)/(1-x**2)
!生成n阶勒让德多项式的导数表达式
End Function g
Real (Kind=DP) Function bis(a, b) !二分法求解函数的解
Implicit None
Real (Kind=DP) :: a, b, c
!a,b是传递进来的划分好的有一个解存在的区间
Do
c = (a+b)/2.0_DP
If (f(c)*f(a)<0) Then
b = c
Else
a = c
End If
If ((b-a)<EPS) exit
End Do
bis = c!bis即是利用二分法求得的解
End Function bis
Subroutine fn0(fn, ak)
Implicit None
Real (Kind=DP) :: m, fn(n), ak(n)
!定义数组,大小n由module开始声明。
Integer :: i, j
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_DP/(n*f1(fn(j))*g(fn(j))) !高斯点的权重
Write (*, *) '高斯点序号', j
Write (*, *) '高斯点', fn(j)
Write (*, *) '高斯点权重', ak(j)
End If
m = m + 0.00001 !执行完一次判断m向前推进一步
End Do
End Subroutine fn0
End Module gsld
Program www_fcode_cn
Use gsld
Implicit None
Real(kind=DP) :: fn(n), ak(n), x, x1, x2, a, b, a1, b1, a2, b2, answer, function_a, function_b, function_c, s1, s2 , p
Integer :: i, j, k
open(55,file="superquadrics_m.inp")
read(55,*) function_a
read(55,*) function_b
read(55,*) function_c
read(55,*) s1
read(55,*) s2
read(55,*) p
close(55)
Call fn0(fn, ak)
a = 0
!积分上下限
a1 = 0
a2 = 0
answer = 0
Do k = 1, n
b2 = function_c
Do j = 1, n
b1 = function_a*(1.0-(x2/function_c)**s1)**(1.0/s1)
Do i = 1, n
b = function_b*(( 1.0-(x2/function_c )**s1)**(s2/s1)-(x1/function_a )**s2)**(1.0/s2)
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))
End Do
End Do
End Do
answer = answer*(b-a)/2.*(b1-a1)/2.*(b2-a2)/2.
Write (*, *) answer
PAUSE
contains
Real(kind=DP) Function y(x, x1, x2) !被积函数
Real(kind=DP) :: x, x1, x2, p
y = 8.* p
End Function y
End Program www_fcode_cn