调用Fortran77中的DDASPK求解器会对数组中的某些点产生奇怪的结果
我用DDASPK用直线法求解抛物型偏微分方程。下面是求解器正在调用的子例程：
SUBROUTINE RES(T,Y,YPRIME,CJ,DELTA,IRES,IPAR,RPAR)
IMPLICIT DOUBLE PRECISION (A-H, O-Z)
DIMENSION Y(*), YPRIME(*), DELTA(*), V(201), LAP3V(201)
H = 0.05D0
NPTS = 201
DO 50 I = 1,NPTS
V(I) = Y(I+1)
50 CONTINUE
DO 60 J = 4,(NPTS-3)
LAP3V(J) = (V(J-3) - 6.0D0*V(J-2) + 15.0D0*V(J-1) -
& 20.0D0*V(J) + 15.0D0*V(J+1) - 6.0D0*V(J+2) +
& V(J+3))/(H**6.0D0)
60 CONTINUE
LAP3V(1) = (2.0D0*V(4) - 12.0D0*V(3) + 30.0D0*V(2) -
&20.0D0*V(1))/(H**6.0D0)
LAP3V(2) = (V(5) - 6.0D0*V(4) + 16.0D0*V(3) - 26.0D0*V(2) +
&15.0D0*V(1))/(H**6.0D0)
LAP3V(3) = (V(6) - 6.0D0*V(5) + 15.0D0*V(4) - 20.0D0*V(3) +
&16.0D0*V(2) - 6.0D0*V(1))/(H**6.0D0)
LAP3V(NPTS-2) = 0.0D0
LAP3V(NPTS-1) = 0.0D0
LAP3V(NPTS) = 0.0D0
DELTA(1) = YPRIME(1) - 1.0D0
DO 70 K = 1,NPTS
DELTA(K+1) = YPRIME(K+1) - LAP3V(K)
70 CONTINUE
RETURN
END这里LAP3V是一维空间中拉普拉斯函数的离散化三次幂(所以是六阶导数)，在第一列中，我们将它与简单的ODE (1)/dt=1耦合在一起(为什么，你会问？当我可以启动并运行它时，我最终将解决更困难的PDE，对于这些PDE，最好将求解器采取的时间步长与解决方案的某些方面耦合起来)。
然而，当我尝试从求解器中调用它时(这是完全隐式的，所以它的第一项工作是在T=0的输入Y的情况下计算一致的YPRIME )，除了两个点之外，一切似乎都正常-在数组的第一个条目中，它很快从0(我在时间0的时候将YPRIME(1) =1作为初始猜测)转换为inf，然后是-nan，原因我不明白-这应该是最简单的部分。同样在YPRIME(NPTS-3)中，我们很快就得到了一个非常大的负数(大约是-1e9)，而它周围的条目是零。我相信这一定与这是在do循环中分配的最终条目有关，而不是逐个案例，但我对fortran还不够熟悉，无法理解修复它的内容或方法。
感谢您的帮助，谢谢。
回答 1
Stack Overflow用户
发布于 2015-03-13 19:56:07
如果没有完整的源代码，就很难确切知道发生了什么。
然而，根据我的经验，每当数组值开始接受奇怪的值(Inf、NaN等)时，几乎总是由于错误地引用了数组值。不正确的数组引用将导致使用引用其他变量的内存，甚至是存储部分二进制可执行文件的位置(在这种情况下，对该位置的赋值通常会导致seg错误和程序崩溃)。
有几种可能的常见场景需要检查(这不一定是详尽的！)：
- 使用假定大小的数组(例如，Y(*))，很容易意外地读取或写入超过数组末尾的数组(例如，在DO 50循环中，如果Y至少有202个元素，那么就有问题了！)。如果你声明了一个可分配的数组，但是分配失败了，请确保你的所有数组都是正确的。
- 你关于你的IPAR和dimensioned.
- If数组被销毁的评论告诉我，其他例程中可能也有一个bug导致了越界的数组赋值。
还要记住，试图求解一组非线性方程来获得一致的导数是非常重要的，并且可能导致找不到有效的解。除了尝试找到更好的初始猜测和/或尝试分析计算导数(如果您可以做到这一点)之外，您对此无能为力。
最后，还有最后一个使用DASPK的技巧。求解器并不关心res子例程的Fortran版本。在Fortran 90中重写RES子例程，从长远来看，您会更快乐！(这也适用于代码的其余部分。Fortran 90比Fortran 77要好得多，除非你别无选择，否则编写新的Fortran 77是疯狂的。) /so 90比Fortran 77好得多。)
https://stackoverflow.com/questions/29022963
复制
腾讯云开发者
