pypinm-koda

# Funkcija Euler
def euler(fun, x, y0, *args, **kwargs):
    """
    Eulerjeva metoda za reševanje navadnih diferencialnih enačb.

    :param fun: Funkcija prvega odvoda y'(x, y)
    :param x: Neodvisna spremenljivka
    :param y0: Začetne vrednosti odvisne spremenljivke y(x=0)
    :return: y: vrednosti odvisne spremenljivke pri podanih vrednosti x
    """
    y = np.zeros((len(x), len(y0)))
    y[0] = np.array(y0).copy()
    for i, xi in enumerate(x[:-1]):
        h = x[i+1] - x[i]
        y[i+1] = y[i] + fun(xi, y[i], *args, **kwargs) * h
    return y.T


# Primer uporabe `scipy.integrate.solve_ivp` z dogodkom:
def upward_cannon(t, y):
    return [y[1], -0.5]
def hit_ground(t, y):
    return y[1]
hit_ground.terminal = True
hit_ground.direction = -1
sol = solve_ivp(upward_cannon, [0, 100], [0, 10], events=hit_ground)
print(sol.t_events)


# Primer reševanje robnega problema (konzola, solve_bvp)
# Podatki:
E = 2.1e5 # N/mm2
I = 5.0e6 # mm^4
l = 1300  # mm
q = 150   # kg/mm

# Funkcija porazdeljene obtežbe:
def obremenitev(x, q=q):
    return np.ones_like(x) * q

# Funkcija prvih odvodov novih spremenljivk - definicija d. e.:
def nosilec(x, y, E=E, I=I, q=q):
    return np.array([y[1], y[2], y[3], obremenitev(x)/E/I])

# Robni preostanki - definicija robnih pogojev:
def robni_preostanki(y_A, y_B):
    return np.array([y_A[0], y_A[1], y_B[0], y_B[2]])

# Neodvisna spremenljivka in začetni približki:
x = np.linspace(0, l, 500)
y0 = np.zeros((4, x.size), dtype=float)

# Rešitev - izračun vpeljanih novih spremeljivk:
rešitev = solve_bvp(nosilec, robni_preostanki, x, y0)