需要对python参数有基本的了解，以便能够编写ODE解决方案的代码系统
需要对python参数有基本的了解，以便能够编写ODE解决方案的代码系统
提问于 2017-07-24 12:14:02
我是一个年长的Fortran程序员，需要帮助一个年轻人使用Heun的方法数值求解一个常微分方程系统。他只知道Python，所以我必须学习最基本的python才能做到这一点。
这是我想出来的。测试代码是一个简单的两自由度系统，每个自由度都有指数增长。
我得到的错误是：
Traceback (most recent call last):
File "program.py", line 28, in <module>
heun (imax, y, dt, t)
File "program.py", line 7, in heun
rhs(y, ydot)
NameError: global name 'ydot' is not defined代码如下：
# Test RHS of ODE system:
def rhs (y, ydot):
ydot[0] = y[0]
ydot[1] = y[1]
return;
# Does one step of Heun's method:
def heun (imax, y, dt, t):
rhs(y, ydot)
for i in range(0, imax):
y_tilde[i] = y[i] + dt * ydot[i]
rhs(y_tilde, ydot_at_tilde)
for i in range(0, imax):
y[i] = y[i] + dt/2 * (ydot[i] + ydot_at_tilde[i])
t = t + dt
return;
# Initial condition
y = [0, 0]
t = 0
dt = 0.01
nsteps = 100
imax = 1
istep = 1
while istep <= nsteps:
heun (imax, y, dt, t)
print istep, y[0], y[1]
istep = istep + 1问:为什么python认为例程heun中的对象ydot是全局的？即使它是全局的，为什么我不能将其作为参数传递？
回答 4
Stack Overflow用户
发布于 2017-07-24 12:26:15
问题出在这里：
def heun(imax, y, dt, t):
rhs(y, ydot)您将使用输入参数y和ydot调用rhs函数。但是在heun函数的作用域中不存在ydot。只有imax、y、dt和t这样做。
类似地，您永远不会定义变量ydot_tilde或ydot_at_tilde
此外，您还需要您的函数返回一些值。
Stack Overflow用户
发布于 2017-07-24 15:36:14
好的，谢谢大家。这是一个代码，它与注释一起工作，以表明我学到了什么：
def rhs (y, ydot):
ydot[0] = y[0]
ydot[1] = y[1]
return ydot;
def heun (ndof, dt, y, t):
# These initializations of local arrays take the place of Fortran declarations:
ydot = [0] * (ndof)
y_tilde = [0] * (ndof)
ydot_at_tilde = [0] * (ndof)
ydot = rhs(y, ydot)
# Note: In python range means:
# range (first element, upto but not including last element)
for i in range(0, ndof):
y_tilde[i] = y[i] + dt * ydot[i]
ydot_at_tilde = rhs(y_tilde, ydot_at_tilde)
for i in range(0, ndof):
y[i] = y[i] + dt/2 * (ydot[i] + ydot_at_tilde[i])
t = t + dt
# Note: This lists the output arguments:
return y, t;
# Initial condition
y = [1, 1]
t = 0
dt = 0.01
nsteps = 100
ndof = 2
istep = 1
while istep <= nsteps:
# Note: This is how you get the output arguments:
y, t = heun (ndof, dt, y, t)
istep = istep + 1
print t, y[0], y[1]Stack Overflow用户
发布于 2017-07-24 15:44:49
这就是我要做的：
import numpy
def heun(ndof, dt, y, t):
ydot = numpy.zeros(ndof)
y_tilde = numpy.zeros(ndof)
ydot_at_tilde = numpy.zeros(ndof)
# Replacing first two elements does not need a function `rhs()`
ydot[:1] = y[:1]
# Vectorized operation, numpy does this loop for you at C speeds
y_tilde = y + dt * ydot
ydot_at_tilde[:1] = y_tilde[:1]
y = y + dt/2 * (ydot + ydot_at_tilde)
t = t + dt
return y, t
y = numpy.ones(2)
t = 0
dt = 0.01
nsteps = 100
ndof = 2
for istep in range(nsteps):
y, t = heun(ndof, dt, y, t)
print(t, y[0], y[1])页面原文内容由Stack Overflow提供。腾讯云小微IT领域专用引擎提供翻译支持
原文链接：
https://stackoverflow.com/questions/45272514
复制相关文章
相似问题
腾讯云开发者
