Как правильно использовать позиционные аргументы для сложной задачи

Question

Я повторно посещаю школьный проект, который я не выполнил к моему удовлетворению. А именно, я написал алгоритм, который принимает ПОЧТИ набор уравнений произвольного размера и решает их итеративно. Проблема в том, что это «почти». По сути, он должен иметь как минимум два уравнения и не будет решать ни одно. Это потому, что я считаю, что я не понимаю, как правильно использовать позиционные аргументы.

Ниже, в основном методе, я определяю две функции 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" в качестве позиционного аргумента, но теперь есть только один, и он больше не повторяется. Итак, я сделал это кортеж 1, когда я передал его в основной метод:

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

Однако, это укусило меня в задницу, потому что теперь я передаю кортеж в мое уравнение 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)
  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

Мое решение в итоге заключалось в том, чтобы преобразовать все списки в пустые массивы, что позволило мне воспользоваться встроенным поэлементным скалярным сложением и умножением. Это сделало вычисление значений «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