[特殊函数] 高斯勒让德求解复杂函数问题

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