如何正确实施scipy.integrate.Radau？

Question

（我什至会感谢 link 的例子，到目前为止我还没有找到。）

我正在尝试使用 scipy.integrate 中的 Radau 来求解二阶微分方程。现在我只是尝试一个简单的例子，这样我就可以理解它是如何工作的（到目前为止不成功）。

假设我有以下等式：

d^2y/dt^2 = C,

这意味着 y = (Ct^2)/2 + Bt.

例如，y(1) = 1 和 C = 2。假设我想找到 t = 10 时的 y 值。

这是我的代码：

from scipy.integrate import Radau
import numpy as np

C = 2.0
y = np.zeros(2)

def fun(t, y):
  #y[0] = C*t
  y[1] = C
  return y

t0 = 1.0
y0 = np.array([1., 1.])
t_bound = 10.0

eq = Radau(fun, t0, y0, t_bound)
print(eq.n)

while(True):
  print(eq.t)
  print(eq.y)
  print(eq.status)
  if eq.status == 'finished':
    break
  eq.step()

输出错误。（如果我取消注释 fun 定义中的注释行，它也会给出错误的答案。但我认为我什至不必告诉求解器，对吧？我通常不知道这个值。）

我认为我最大的问题是我不确定应该将什么传递为 fun。文档说它应该在系统的右侧，所以我认为一阶推导应该在 y[0]，二阶推导在 y[1] 等

我做错了什么？应该如何实施？

Answer 1

Atm 你正在解决

y0' = y0,  y0(1)=1
y1' = 2,   y1(1)=1

有解决方案 y0(t)=exp(t-1) 和 y1(t)=2*t-1 这肯定不是你想要的。你想要一阶系统

y0' = y1
y1' = C

所以你需要

def fun(t,y): return [y[1], C]

则解为y1(t)=C*t+B=2*t-1和y0(t)=0.5*C*t^2+B*t+A=t^2-t+1，积分正确结束为eq.y = [91. 19.]。

如何正确实施scipy.integrate.Radau？

How to properly implement scipy.integrate.Radau?

python

scipy

differential-equations