1

I was trying to solve a system of 2 differential equations using python script, but it seems I'm doing something wrong.

I have 2 equations:

dx/dt = ax - bxy
dy/dt = -cx + dxy

The integration interval is [15; 45]. a, b, c, d coefficients are given as well.

The code I`m using to calculate this is below:

def get_next_v(_v_i: float, _q: List[float]) -> float:
    assert len(_q) == 4
    _q_s = _q[0] + 2 * _q[1] + 2 * _q[2] + _q[3]
    return _v_i + ((1 / 6) * _q_s)

def get_qi_seq(_h: float, _t_i: float, _x_i: float, _y_i: float,
               _a: float, _b: float, _c: float, _d: float) -> Tuple[List[float], List[float]]:

    # x' = ax - bxy
    # y` = -cy + dxy

    _k_1 = _h * x(_x_i, _y_i, _a, _b)
    _q_1 = _h * y(_x_i, _y_i, _c, _d)

    _k_2 = _h * x(_x_i + _h / 2, _y_i + _k_1 / 2, _a, _b)
    _q_2 = _h * y(_x_i + _h / 2, _y_i + _q_1 / 2, _c, _d)

    _k_3 = _h * x(_x_i + _h / 2, _y_i + _k_2 / 2, _a, _b)
    _q_3 = _h * y(_x_i + _h / 2, _y_i + _q_2 / 2, _c, _d)

    _k_4 = _h * x(_x_i + _h, _y_i + _k_3, _a, _b)
    _q_4 = _h * y(_x_i + _h, _y_i + _q_3, _c, _d)

    return [_k_1, _k_2, _k_3, _k_4], [_q_1, _q_2, _q_3, _q_4]

def x(_x_i: float, _y_i: float, _a: float, _b: float) -> float:
    return (_x_i * _a) - (_b * _x_i * _y_i)


def y(_x_i: float, _y_i: float, _c: float, _d: float) -> float:
    return (- (_y_i * _c)) + (_d * _x_i * _y_i)


def get_x_y_vals(_t0: float, _T: float, _x0: float, _y0: float, _h: float,
                 _a: float, _b: float, _c: float, _d: float) -> List[List[float]]:

    _t = _t0
    _xi = _x0
    _yi = _y0

    _res = [[_t, _xi, _yi]]
    while _t < _T:

        _ki, _qi = get_qi_seq(_h, _t, _xi, _yi, _a, _b, _c, _d)

        _xi = get_next_v(_xi, _ki)
        _yi = get_next_v(_xi, _qi)
        _t += _h
        _res.append([_t, _xi, _yi])

    return _res

def main():
    # # a = 1.75
    # # b = 2.15
    # # c = 1.35
    # # d = 0.82
    # # interval = [14, 44]
    # # e = 10 ** -3
    #
    a = 1.77
    b = 2.17
    c = 1.38
    d = 0.89
    x0 = 3.39
    y0 = 2.13
    interval = [15, 45]
    e = 10 ** -3

    vals = get_x_y_vals(interval[0],
                        interval[1],
                        x0,
                        y0,
                        0.1,
                        a, b, c, d)

    print(get_h(e))
    print(vals)

    t_range = [_item[0] for _item in vals]
    x_range = [_item[1] for _item in vals]
    y_range = [_item[2] for _item in vals]

    plt.plot(t_range, x_range)
    plt.plot(t_range, y_range)
    # plt.plot(x_range, y_range)
    plt.show()

The problem is, that the example from the book I`m trying to replicate has an absolutely different plot.

Example plot:

Example plot

My result:

My plot

What am I doing wrong and how can I fix this?

MattDMo
  • 100,794
  • 21
  • 241
  • 231
Anton Grant
  • 23
  • 7
  • The updates for the `_x` are always the `_k` named result of the `x()` function, for `y` accordingly `_q`. There should be no mixing of these associations. See https://stackoverflow.com/questions/58255653/lotka-volterra-with-runge-kutta-not-desired-precision for a similar (but not equal) problem. – Lutz Lehmann Dec 18 '22 at 09:55
  • 1
    @LutzLehmann thank you for your help, the answer you linked helped a lot! – Anton Grant Dec 18 '22 at 12:20

0 Answers0