如何为复杂问题恰当地使用位置参数

Question

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

下面，在main方法中，我定义了两个函数y\u prime和z\u 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”作为一个位置参数，但现在只有一个，它不再是iterable。因此，当我在main方法中传递它时，我将它设为1的元组：

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

不过，这对我来说有点不对劲，因为现在我把一个元组传递到我的y\u素数方程中，而它真正需要的是一个浮点数：

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