高斯勒让德求解复杂函数问题

无奈的小强

[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 

kerb

http://bbs.fcode.cn/forum.php?mo ... &extra=#pid5135
建议参考这个帖子，
如果希望精度更高，可以选用更高阶数，或者增加区间分段数
一般情况下，是对所给区间分段，每一段分别使用高斯--勒让德积分，最后加起来，这样精度才会高，要避免在所给整个区间直接使用高斯--勒让德积分

pasuka

kerb 发表于 2016-8-17 21:07
http://bbs.fcode.cn/forum.php?mod=viewthread&tid=814&page=1&extra=#pid5135
建议参考这个帖子，
如果 ...

lz利用多项式求根方法，是否比较过60阶之后的结果？印象中十多年前有老外在mathworks网站上贴出对比，多项式求根方法阶次过60后，积分权重出现负值。一般而言，低于10阶，查手册，高的话构造伴随矩阵算特征值，类似的matlab代码很多，改回fortran不难。*************************************************************************************************************************
补充一点干货
ACM的CALGO下面的Fortran现成代码
http://people.sc.fsu.edu/~jburkardt/f_src/toms655/toms655.html
某老外几年前的代码且附带高斯积分点结果
http://www.holoborodko.com/pavel ... erical-integration/

无奈的小强

kerb 发表于 2016-8-17 21:07
http://bbs.fcode.cn/forum.php?mod=viewthread&tid=814&page=1&extra=#pid5135
建议参考这个帖子，
如果 ...

您好，我想问下，我在fortran中用拟牛顿broyden法求解非线性超函数的方程组时，总是出现NAN情况，您知道是为什么吗？

无奈的小强

pasuka 发表于 2016-8-17 21:41
lz利用多项式求根方法，是否比较过60阶之后的结果？印象中十多年前有老外在mathworks网站上贴出对比，多项 ...

您好，我想问下，我在fortran中用拟牛顿broyden法求解非线性超函数的方程组时，总是出现NAN情况，您知道是为什么吗？

kerb

pasuka 发表于 2016-8-17 21:41
lz利用多项式求根方法，是否比较过60阶之后的结果？印象中十多年前有老外在mathworks网站上贴出对比，多项 ...

我不明白你的意思

pasuka

kerb 发表于 2016-8-23 22:05
我不明白你的意思

多项式直接二分法求根的办法，稳定性随着阶次增加会下降
构造伴随矩阵求解能够避免上述问题，这个是计算方法课程的基本内容。。。

kerb

pasuka 发表于 2016-8-24 08:28
多项式直接二分法求根的办法，稳定性随着阶次增加会下降
构造伴随矩阵求解能够避免上述问题，这个是计算 ...

我不想与你讨论这些问题，请打住！

pasuka

算法研究版面不讨论算法的稳定性、收敛性、鲁棒性等等有趣的问题，那就不好玩了。。。

kerb

pasuka 发表于 2016-8-24 12:29
算法研究版面不讨论算法的稳定性、收敛性、鲁棒性等等有趣的问题，那就不好玩了。。。 ...

说句刺激一点的话，叫做自以为是！