ODE 常微分方程组求解
大概的一个思路
- 首先考虑一阶方程: $\dot{y} = f(y, t)$, 构造迭代器, 如 $y_{n+1} = y_n + f(y_n, t_n) \mathrm{d} t$ (Forward Euler), 根据初始值迭代即可得到一条路径
- 考虑多元一阶方程, 相当于是将 1 中的 $y$ 和 $f$ 拓展为一个向量
ndarray
- 考虑二阶方程: $\ddot{y} = f(y, \dot{y}, t)$, 引入变元 $u = \dot{y}, \ddot{y} = u \dot{u}$, 添加关于 $u$ 的初始条件, 故问题退化为 2 的多元一阶方程
Runge-Kutta
这里提供的是 Explicit Runge-Kutta 的方法, 想要了解 Forward Euler 或者 Backward Euler, 或者是 Implicit Runge-Kutta 的方法, 请参考 Runge-Kutta(Wikipedia).
Runge-Kutta 的一般形式如下:
$$y_{n+1} = y_n + h \sum b_i k_i$$
其中: $k_i = f(t_n + c_i h, y_n + (a{i,1} k_1 + \dots + a{i,i-1} k_{s-1}))$.
这里是一个简单的解释: 考虑之前的 数值求导的多级展开, 可以简单理解为这里的 $k_i$ 就是多级展开中的展开项, 是为提供 $\dot{y}$ 更好的近似.
于是可以构造一个抽象的迭代器:
def make_runge_kutta_iter(bs, c_as):
def iter(f, y, time, dt):
route = [y]
t = 0
while t < time:
ks = []
for (c, *a_s) in c_as:
sum_ks = sum([a * k for (a, k) in zip(a_s, ks)])
ks.append(f(t + c * dt, y + sum_ks * dt))
y = y + dt * sum([b * k for (b, k) in zip(bs, ks)])
route.append(y)
t += dt
return route
return iter
注: 这里我是从 Lisp 代码修改过来的, 还来不及 debug, 可能存在一些小问题, 我会在最后放上我的 Lisp 代码 (这是完全可以跑的)
[来自下午的更新]: 修改了输入的形式, 让构造过程更加的简单一点; 发现一开始的表格输错了, 尴尬, 还以为是代码不对.
于是考虑 Forward Euler, 也就是 1st Runge-Kutta:
如下所示:
rk1 = make_runge_kutta_iter([1],
[[0]])
或者 4th Runge-Kutta:
c | as | as | as | as |
0 | | | | |
1/2 | 1/2 | | | |
1/2 | 0 | 1/2 | | |
1 | 0 | 0 | 1 | |
----- | ----- | ----- | ----- | ----- |
| 1/6 | 1/3 | 1/3 | 1/6 |
rk4 = make_runge_kutta_iter([1/6, 1/3, 1/3, 1/6],
[[0],
[1/2, 1/2],
[1/2, 0, 1/2],
[1, 0, 0, 1]])
理论上来说适合任意的 Runge Kutta 函数, 这里再多抄几个 Wikipedia 上的结果, 在 Lisp 上测试的结果是看着都差不多.
rk4_2 = make_runge_kutta_iter([1/8, 3/8, 3/8, 1/8],
[[0],
[1/3, 1/3],
[2/3, -1/3, 1]
[1, 1, -1, 1]])
运行迭代器:
def eqs(t, arg, k = 1.0, tau = 1.0):
x, v = arg
return np.array([
v, # x' = v
- k * x - tau * v # v' = - k * x - tau * v
])
rk4(eqs, np.array([1, 0]), 10, 0.1)
关于方法的选择和 dt
的选择
首先, 一般来说不会选择 Forward Euler 或者 Backward Euler, 因为效果不佳.
然后, 随着迭代的时间 time
的增加, 误差一定会越来越大, 想要减少误差, 就要减少 dt
, 一般来说, 减少 dt
比如 10 倍, 那么按照 $O(hn)$, 对应可以增加的迭代时间也可以增加 $10n$ 倍.
Lisp Code
下面是我的 Lisp 的代码:
(defun make-runge-kutta-iter (bs c-as)
(let ((bs (mapcar #'float bs)) ; trun to float type for fast computing
(c-as (mapcar (lambda (c-a) (mapcar #'float c-a)) c-as)))
(macrolet ((w-sum (ws ks)
`(reduce #'add (mapcar #'mul ,ws ,ks)
:initial-value 0)))
(lambda (fn init-y times dt)
(loop for y = init-y ; y + dt * sum(bi * ki)
then (add y (mul dt (w-sum bs ks)))
for time below times by dt
;; ki = f(t + ci * dt, y + sum(a_ij * k_j) * dt)
for ks = (loop for (c . as) in c-as
collect (apply fn (+ time (* c dt))
(add y (mul (w-sum as ks) dt)))
into ks
finally (return ks))
collect y)))))
(defmacro def-runge-kutta (name bs &body c-as)
"Define a explicit Runge-Kutta process marco `name' and given tableau.
The `bs' will be a list like (b1 b2 ... bs);
The `c-as' will be like (c a1 a2 ... as-1).
Example:
(def-runge-kutta rk1
(1)
(0))"
(with-gensyms (iter)
`(let ((,iter (make-runge-kutta-iter ',bs ',c-as)))
(defmacro ,name (lambda-list fn-body
&key (dt 0.001) (time 1.0) (var-t 'ti))
(let* ((ys (mapcar #'first lambda-list))
(y0 (mapcar #'second lambda-list))
(body (if (eq (car fn-body) 'quote)
(second fn-body)
(cons 'list fn-body))))
`(funcall ,,iter (lambda (,var-t ,@ys)
(declare (ignorable ,var-t))
,body)
(list ,@y0) ,time ,dt))))))
;; 1st Runge-Kutta
(def-runge-kutta rk1 (1) (0))
;; 2nd Runge-Kutta
(def-runge-kutta rk2
(0 1)
(0)
(1/2 1/2))
;; 4th Runge-Kutta
(def-runge-kutta rk4
(1/6 1/3 1/3 1/6)
(0)
(1/2 1/2)
(1/2 0 1/2)
(1 0 0 1))
其中用到的一些依赖代码:
(macrolet ((listfy-fn (name fn)
`(defun ,name (a b)
(cond ((and (numberp a) (numberp b))
(,fn a b))
((and (numberp a) (listp b))
(mapcar (lambda (v) (,name a v)) b))
((and (listp a) (numberp b))
(mapcar (lambda (v) (,name v b)) a))
(t (mapcar #',name a b))))))
(listfy-fn add +)
(listfy-fn sub -)
(listfy-fn mul *)
(listfy-fn div /))
(defun sum (&rest vars)
(reduce #'add vars :initial-value 0))