高斯勒让德求解复杂函数问题
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
http://bbs.fcode.cn/forum.php?mod=viewthread&tid=814&page=1&extra=#pid5135
建议参考这个帖子,
如果希望精度更高,可以选用更高阶数,或者增加区间分段数
一般情况下,是对所给区间分段,每一段分别使用高斯--勒让德积分,最后加起来,这样精度才会高,要避免在所给整个区间直接使用高斯--勒让德积分 本帖最后由 pasuka 于 2016-8-22 08:50 编辑
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/numerical-methods/numerical-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情况,您知道是为什么吗? pasuka 发表于 2016-8-17 21:41
lz利用多项式求根方法,是否比较过60阶之后的结果?印象中十多年前有老外在mathworks网站上贴出对比,多项 ...
我不明白你的意思 kerb 发表于 2016-8-23 22:05
我不明白你的意思
多项式直接二分法求根的办法,稳定性随着阶次增加会下降
构造伴随矩阵求解能够避免上述问题,这个是计算方法课程的基本内容。。。 本帖最后由 kerb 于 2016-8-24 12:09 编辑
pasuka 发表于 2016-8-24 08:28
多项式直接二分法求根的办法,稳定性随着阶次增加会下降
构造伴随矩阵求解能够避免上述问题,这个是计算 ...
我不想与你讨论这些问题,请打住! 算法研究版面不讨论算法的稳定性、收敛性、鲁棒性等等有趣的问题,那就不好玩了。。。 pasuka 发表于 2016-8-24 12:29
算法研究版面不讨论算法的稳定性、收敛性、鲁棒性等等有趣的问题,那就不好玩了。。。 ...
说句刺激一点的话,叫做自以为是!
页:
[1]
2