下面的代码使用正割方法来查找解析函数的根。解析函数f必须在代码的函数部分指定。下面的代码运行良好，没有编译错误。但是，对于我要解决的问题，我不知道解析函数f。
相反，我用数字计算函数，并将其存储为数组。我想现在应用我的代码来找到这个函数的根。所以，我如何修改我的代码，使输入不是一个解析函数，而只是一个数组，我已经计算？
我的工作代码如下，我想我只需要修改最后一部分，我调用函数f，我只是不确定如何去做。谢谢！
program main
implicit none
real :: a = 1.0, b = -1.0
integer :: m = 8
interface
function f(x)
real, intent(in) :: x
end function
end interface
call secant(f,a,b,m)
end program main
subroutine secant(f,a,b,m)
implicit none
real, intent(in out) :: a,b
integer, intent(in) :: m
real :: fa, fb, temp
integer :: n
interface
function f(x)
real, intent(in) :: x
end function f
end interface
fa = f(a)
fb = f(b)
if (abs(fa) > abs(fb)) then
temp = a
a = b
b = temp
temp = fa
fa = fb
fb = temp
end if
print *," n x(n) f(x(n))"
print *," 0 ", a, fa
print *," 1 ", b, fb
do n = 2,m
if (abs(fa) > abs(fb)) then
temp = a
a = b
b = temp
temp = fa
fa = fb
fb = temp
end if
temp = (b - a)/(fb - fa)
b = a
fb = fa
a = a - fa*temp
fa = f(a)
print *,n,a,fa
end do
end subroutine secant
real function f(x)
implicit none
real, intent(in) :: x
f = x**5 + x**3 + 3.0 !analytic form of a function, I don't actually have this though, I just have the function stored as an array
end function f
最佳答案
我想在我的评论中说的如下。
您可以修改您的secant子例程以获取一个抽象类（FAZ）的对象，该类保证有一个函数f。例如，如下所示。solver.f90
!*****************************************************************
MODULE solver
!*****************************************************************
IMPLICIT NONE
PRIVATE
PUBLIC FAZ
PUBLIC secant
TYPE, ABSTRACT :: FAZ
CONTAINS
PROCEDURE(f), deferred, pass :: f
END TYPE FAZ
ABSTRACT INTERFACE
FUNCTION f(this, x)
IMPORT :: FAZ
REAL :: f
CLASS(FAZ), INTENT(IN) :: this
REAL, INTENT(IN) :: x
END FUNCTION f
END INTERFACE
!=====================================================================
CONTAINS
!=====================================================================
subroutine secant(oFAZ,a,b,m)
CLASS(FAZ) :: oFAZ
real, intent(in out) :: a,b
integer, intent(in) :: m
real :: fa, fb, temp
integer :: n
fa = oFAZ%f(a)
fb = oFAZ%f(b)
if (abs(fa) > abs(fb)) then
temp = a
a = b
b = temp
temp = fa
fa = fb
fb = temp
end if
print *," n x(n) f(x(n))"
print *," 0 ", a, fa
print *," 1 ", b, fb
do n = 2,m
if (abs(fa) > abs(fb)) then
temp = a
a = b
b = temp
temp = fa
fa = fb
fb = temp
end if
temp = (b - a)/(fb - fa)
b = a
fb = fa
a = a - fa*temp
fa = oFAZ%f(a)
print *,n,a,fa
end do
end subroutine secant
END MODULE solver
然后，通过将抽象类
f扩展到具体类FAZ，可以以任何方式实现函数MyFAZ的行为。例如，我写的如下。myfaz.f90!*******************************************************************
MODULE my_concrete_faz
!*******************************************************************
USE solver, ONLY : FAZ
IMPLICIT NONE
PRIVATE
PUBLIC MyFAZ
PUBLIC MyFAZ_constructor
TYPE, EXTENDS(FAZ) :: MyFAZ
PRIVATE
REAL, DIMENSION(:), ALLOCATABLE :: xdata, fdata
CONTAINS
PROCEDURE :: destructor
PROCEDURE :: f
END TYPE MyFAZ
! ================================================================
CONTAINS
! ================================================================
! ****************************************************************
FUNCTION MyFAZ_constructor(xdata_arg, fdata_arg) RESULT(oMyFAZ)
! ****************************************************************
TYPE(MyFAZ) :: oMyFAZ
REAL, DIMENSION(:), INTENT(IN) :: xdata_arg, fdata_arg
INTEGER :: ndata, jj
ndata = size(xdata_arg)
if (size(fdata_arg) /= ndata) then
stop 'MyFAZ_constructor: array size mismatch .. ndata'
end if
do jj=1,ndata-1
if (xdata_arg(jj)>xdata_arg(jj+1)) then
stop 'MyFAZ_constructor: expecting a sorted xdata. I am lazy.'
end if
end do
allocate(oMyFAZ%xdata(ndata))
allocate(oMyFAZ%fdata(ndata))
oMyFAZ%xdata = xdata_arg
oMyFAZ%fdata = fdata_arg
END FUNCTION MyFAZ_constructor
! ****************************************************************
SUBROUTINE destructor(this)
! ****************************************************************
CLASS(MyFAZ), INTENT(INOUT) :: this
deallocate(this%xdata)
deallocate(this%fdata)
END SUBROUTINE destructor
! ****************************************************************
FUNCTION f(this, x)
! ****************************************************************
! evaluates the function.
! Linear interpolation is used here, but this will not make sense
! in actual application. Everything is written in a very inefficient way.
REAL :: f
CLASS(MyFAZ), INTENT(IN) :: this
REAL, INTENT(IN) :: x
!
INTEGER :: jj
REAL :: rr
do jj=1, size(this%xdata)-1
if (this%xdata(jj)<=x .and. x<=this%xdata(jj+1)) then
exit
end if
end do
rr = (this%fdata(jj+1) - this%fdata(jj))/(this%xdata(jj+1) - this%xdata(jj))
f = rr*(x - this%xdata(jj)) + this%fdata(jj)
END FUNCTION f
END MODULE my_concrete_faz
我用线性插值，只是为了演示。实际上，如果
f(x) = r x + s，那么不用正割方法就知道解。您将有自己合适的方法来评估数据点之间的
f(x)。您可以使用以上两个模块，如下所示。
main.f90PROGRAM demo
USE solver, ONLY : secant
USE my_concrete_faz, ONLY : MyFAZ, MyFAZ_constructor
IMPLICIT NONE
REAL, DIMENSION(:), ALLOCATABLE :: xdata, fdata
INTEGER :: ndata
INTEGER :: niter_max
REAL :: xa, xb
TYPE(MyFAZ) :: oMyFAZ
niter_max = 10
xa = -2.0
xb = 3.0
! prepare data
ndata = 4
allocate(xdata(ndata))
allocate(fdata(ndata))
xdata(1) = -3.0
xdata(2) = -1.1
xdata(3) = 1.2
xdata(4) = 3.8
fdata(1) = -1.5
fdata(2) = -0.9
fdata(3) = 0.1
fdata(4) = 0.8
! prepare the function
oMyFAZ = MyFAZ_constructor(xdata, fdata)
deallocate(xdata)
deallocate(fdata)
! solve
call secant(oMyFAZ,xa,xb,niter_max)
write(*,*) '**************'
write(*,*) 'normal end'
write(*,*) '**************'
END PROGRAM demo
我编译、构建并得到如下输出。
$ ifort -c solver.f90
$ ifort -c myfaz.f90
$ ifort -c main.f90
$ ifort -o demo *.o
$ ./demo
n x(n) f(x(n))
0 3.000000 0.5846154
1 -2.000000 -1.184211
2 1.347448 0.1396975
3 0.8285716 -6.1490655E-02
4 0.9871597 7.4606538E-03
5 0.9700001 0.0000000E+00
6 0.9700001 0.0000000E+00
7 NaN NaN
8 NaN NaN
9 NaN NaN
10 NaN NaN
**************
normal end
**************
$
因为您的
NaN子程序在最大迭代之前到达了解决方案，但是在循环的中间没有出路。这是一张数据图。