- html - 出于某种原因,IE8 对我的 Sass 文件中继承的 html5 CSS 不友好?
- JMeter 在响应断言中使用 span 标签的问题
- html - 在 :hover and :active? 上具有不同效果的 CSS 动画
- html - 相对于居中的 html 内容固定的 CSS 重复背景?
我正在使用 SciPy 的 integrate.odeint 函数集成刚性 ODE 系统。由于集成非常重要且耗时,我还使用了相应的 jacobian。通过重新排列方程,我可以将雅可比矩阵定义为带状矩阵。关注API documentation我想用 mu 和 ml 参数定义形状。不幸的是,文档有点含糊不清,因此我无法弄清楚如何实现我的 jacobian 函数。
为了验证必须如何调用 odeint,我一直在使用以下(有点傻)代码:
from scipy.integrate import odeint
lmax = 5
def f1(y, t):
ydot = np.zeros(lmax)
for i in range(lmax):
ydot[i] = y[i]**2-y[i]**3
return ydot
def fulljac(y, t,):
J = np.zeros((lmax, lmax))
J[0,0]=-3*y[0]**2 + 2*y[0]
J[1,1]=-3*y[1]**2 + 2*y[1]
J[2,2]=-3*y[2]**2 + 2*y[2]
J[3,3]=-3*y[3]**2 + 2*y[3]
J[4,4]=-3*y[4]**2 + 2*y[4]
return J
## initial conditions and output times
delta = 0.0001;
yini = np.array([delta]*lmax)
times = np.linspace(0, 2/delta, 100)
y, infodict = odeint(f1, yini, times, Dfun=fulljac, full_output=1)
print("f1: nst: {0}, nfe: {1}, nje: {2}".format(infodict["nst"][-1],
infodict["nfe"][-1],
infodict["nje"][-1]))
使用完整的 NxN 雅可比矩阵集成成功。仅使用对角线和 mu=0 和 ml=0 积分也成功。
为了测试带状矩阵用例,我使用 mu=1 和 ml=1 创建了一个人工 3xN 带状雅可比矩阵,其中对角线外的所有导数都是零。这会导致求解器出现奇怪的行为(类似于我在原始问题中看到的非对角线不为零的情况)。
def bandjac(y, t):
J = np.zeros((lmax, 3))
J[0,1]=-3*y[0]**2 + 2*y[0]
J[1,1]=-3*y[1]**2 + 2*y[1]
J[2,1]=-3*y[2]**2 + 2*y[2]
J[3,1]=-3*y[3]**2 + 2*y[3]
J[4,1]=-3*y[4]**2 + 2*y[4]
return J
y, infodict = odeint(f1, yini, times, Dfun=bandjac, full_output=1, mu=1, ml=1)
print("f1: nst: {0}, nfe: {1}, nje: {2}".format(infodict["nst"][-1],
infodict["nfe"][-1],
infodict["nje"][-1]))
在 SciPy 的 integrate.odeint 中使用 banded jacobian 选项的正确方法是什么?
最佳答案
为了完整起见,我正在回答我自己的问题。
正如 Warren Weckesser 指出的那样,有一个 bug在 Scipy <0.14.0 中关于 odeint 如何处理带状雅可比矩阵。
odeint 状态的当前文档:
Dfun should return a matrix whose rows contain the non-zero bands (starting with the lowest diagonal).
我认为这是不正确的。相反,它应该从最高对角线开始。
以下代码片段显示了 Dfun 应如何返回雅可比矩阵(源自 test_integrate.py unit test ):
def func(y, t, c):
return c.dot(y)
def jac(y, t, c):
return c
def bjac_rows(y, t, c):
return np.array([[ 0, 75, 1, 0.2], # mu - upper
[ -50, -0.1, -1e-04, -40], # diagonal
[ 0.2, 0.1, -0.2, 0]]) # lower - ml
c = array([[-50, 75, 0, 0],
[0.2, -0.1, 1, 0],
[0, 0.1, -1e-4, 0.2],
[0, 0, -0.2, -40]])
y0 = arange(4)
t = np.linspace(0, 50, 6)
# using the full jacobian
sol0, info0 = odeint(func, y0, t, args=(c,), full_output=True,Dfun=jac)
print("f1: nst: {0}, nfe: {1}, nje: {2}".format(info0["nst"][-1],
info0["nfe"][-1],
info0["nje"][-1]))
# using the row based banded jacobian
sol2, info2 = odeint(func, y0, t, args=(c,), full_output=True,Dfun=bjac_rows, ml=1, mu=1)
print("f1: nst: {0}, nfe: {1}, nje: {2}".format(info2["nst"][-1],
info2["nfe"][-1],
info2["nje"][-1]))
注意:转置带状矩阵似乎不适用于col_deriv=True
关于python - 带有带状雅可比矩阵的 Scipy odeint,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/25978099/
我对 Elm 进行排序时使用什么排序算法 List ? > sort [1,3,5,6] [1,3,5,6] : [comparable] 什么是 [comparable] 类型以及如何取回 numb
我必须编写一个优先队列作为以下接口(interface)的实现: public interface PQueue> { public void insert( T o ); // insert
设以下实体: @Entity public class Person { @Id long id; @ManyToOne Family fam; @ManyTo
今天在 AP 计算机科学课上,我有这段代码: Comparable x = 45; Comparable y = 56; System.out.println(x.compar
如果您知道 WPF 的 MVVM 模式,那么您就会知道 Josh smith msdn 文章,其中 CustomerViewModel 不包含如下简单属性: public string FirstNa
我是一名优秀的程序员,十分优秀!