python - 如何针对复杂问题适本地使用位置参数

Question

我正在重新审视一个学校项目，但我对它的完成并不满意。也就是说，我编写了一个算法，它采用几乎任意大小的方程组并迭代求解它们。问题是“几乎”的部分。本质上，它必须至少有两个方程，并且不会解出一个方程。这是因为，我相信，我不明白如何正确使用位置参数。

下面，在 main 方法中，我定义了两个函数 y_prime 和 z_prime。如果我都通过了它们，我会得到一张漂亮的解决方案图表。但是，如果我只将 y_prime 及其初始条件和解向量传递给 rungekutta() 函数，事情就会变得一团糟:

import numpy as np
import matplotlib.pyplot as plt
import pandas as pd


def rungekutta(dt, y, t, *funcs):
    """
    The following code was written in order to
    reproduce the classic 4th order Runge-Kutta numerical
    method of solving a system of differential equations.
    The aim was to not only apply this to the budworm deforestation
    model developed by Ludwig et al, but also to create an algorithm
    that is generic enough to accept a wide range of ODEs and
    systems of ODEs.

    :param dt: time step "Delta t"
    :param y: The solution vector at the last time step
    :param t: The time at the last time step
    :param funcs: the vector field dy/dt = f(t,y)
    :return: The solution vector for the next time step
    """

    k1 = [dt * f(*y, t) for f in funcs]
    args = [y_n + 0.5 * k_1 for y_n, k_1 in zip((*y, t), (*k1, dt))]
    k2 = [dt * f(*args) for f in funcs]
    args = [y_n + 0.5 * k_2 for y_n, k_2 in zip((*y, t), (*k2, dt))]
    k3 = [dt * f(*args) for f in funcs]
    args = [y_n + k_3 for y_n, k_3 in zip((*y, t), (*k3, dt))]
    k4 = [dt * f(*args) for f in funcs]

    return [y_n + (k_1 + 2 * k_2 + 2 * k_3 + k_4) / 6 for y_n, k_1, k_2, k_3, k_4 in
            zip(y, k1, k2, k3, k4)]


if __name__ == '__main__':

    def y_prime(y, z, t):
        return -t * y

    def z_prime(y, z, t):
        return z

    t_0 = -10
    t_n = 10
    dt = .05

    steps = int((t_n - t_0) / dt)

    y_soln = [0] * steps
    z_soln = [0] * steps
    time = np.arange(t_0, t_n, dt)

    y_soln[0] = 1.928749848e-22
    z_soln[0] = .0000453999297625

    for i in np.arange(1, steps):
        y_soln[i] = rungekutta(dt, y_soln[i-1], time[i-1], y_prime)

我在尝试传递单个方程式时收到的第一个错误是:

Traceback (most recent call last):
  File "C:/Users/wesle/PycharmProjects/Budworms/RK4v2.py", line 57, in <module>
    y_soln[i] = rungekutta(dt, y_soln[i-1], time[i-1], y_prime, z_prime)
  File "C:/Users/wesle/PycharmProjects/Budworms/RK4v2.py", line 23, in rungekutta
    k1 = [dt * f(*y, t) for f in funcs]
  File "C:/Users/wesle/PycharmProjects/Budworms/RK4v2.py", line 23, in <listcomp>
    k1 = [dt * f(*y, t) for f in funcs]
TypeError: y_prime() argument after * must be an iterable, not float

这是因为，我想，我有“y_soln”作为位置参数，但现在只有一个，而且它不再是可迭代的。因此，当我在 main 方法中传递它时，我将它设为 1 的元组:

for i in np.arange(1, steps):
    y_soln[i] = rungekutta(dt, (y_soln[i-1],), time[i-1], y_prime)

然而，这让我很生气，因为现在我将一个元组传递到我的 y_prime 方程中，而它真正需要的是一个 float :

Traceback (most recent call last):
  File "C:/Users/wesle/PycharmProjects/Budworms/RK4v2.py", line 57, in <module>
    y_soln[i] = rungekutta(dt, (y_soln[i-1],), time[i-1], y_prime)
  File "C:/Users/wesle/PycharmProjects/Budworms/RK4v2.py", line 23, in rungekutta
    k1 = [dt * f(*y, t) for f in funcs]
  File "C:/Users/wesle/PycharmProjects/Budworms/RK4v2.py", line 23, in <listcomp>
    k1 = [dt * f(*y, t) for f in funcs]
  File "C:/Users/wesle/PycharmProjects/Budworms/RK4v2.py", line 38, in y_prime
    return -t * y
TypeError: can't multiply sequence by non-int of type 'numpy.float64'

到目前为止，我唯一的变通方法是除了求解我感兴趣的方程之外，还求解一个额外的随机方程，例如 $y= y'$。尽管这看起来效率很低。

所以，似乎我这样做是该死的，如果我不这样做也是该死的。有什么补救办法吗？

编辑如果你想看到代码实际工作，替换这个:

  for i in np.arange(1, steps):
        y_soln[i] = rungekutta(dt, (y_soln[i-1],), time[i-1], y_prime)

在我将方程及其解向量传递给函数的实例中:

for i in np.arange(1, steps):
    y_soln[i], z_soln[i] = rungekutta(dt, (y_soln[i-1], z_soln[i-1]), time[i-1], y_prime, z_prime)

Answer 1

我的解决方案最终是将所有列表转换为 numpy 数组，这使我能够利用内置的逐元素标量加法和乘法。这使得计算“k”值变得不那么麻烦和复杂:

def rk4(dt, t, field, y_0):
    """
    :param dt: float - the timestep
    :param t: array - the time mesh
    :param field: method - the vector field y' = f(t, y)
    :param y_0: array - contains initial conditions
    :return: ndarray - solution
    """

    # Initialize solution matrix. Each row is the solution to the system
    # for a given time step. Each column is the full solution for a single
    # equation.
    y = np.asarray(len(t) * [y_0])

    for i in np.arange(len(t) - 1):
        k1 = dt * field(t[i], y[i])
        k2 = dt * field(t[i] + 0.5 * dt, y[i] + 0.5 * k1)
        k3 = dt * field(t[i] + 0.5 * dt, y[i] + 0.5 * k2)
        k4 = dt * field(t[i] + dt, y[i] + k3)
        y[i + 1] = y[i] + (k1 + 2 * k2 + 2 * k3 + k4) / 6

    return y