Решение нелинейных дифференциальных уравнений первого порядка с использованием Python

Question

Я хотел бы решить нелинейное дифференциальное уравнение первого порядка с использованием Python.

Например,

ДФ/дт = е ** 4

Я написал следующую программу, но у меня есть проблема с Matplotlib, так что я не знаю, если метод, который я использовал с SciPy верно.

from scipy.integrate import odeint 
import numpy as np 
import matplotlib.pyplot as plt 
derivate=lambda f,t: f**4 
f0=10 
t=np.linspace(0,2,100) 
f_numeric=scipy.integrate.odeint(derivate,f0,t) 
print(f_numeric) 
plt.plot(t,f_numeric) 
plt.show() 

что приводит к следующей ошибке:

AttributeError: 'float' object has no attribute 'rint' 

Answer 1

В этом случае, вы могли бы быть лучше использовать Sympy, что позволяет получить замкнутую форму решения:

from IPython.display import display 
import sympy as sy 
from sympy.solvers.ode import dsolve 
import matplotlib.pyplot as plt 
import numpy as np 

sy.init_printing() # LaTeX like pretty printing for IPython 


t = sy.symbols("t", real=True) 
f = sy.symbols("f", function=True) 


eq1 = sy.Eq(f(t).diff(t), f(t)**4) # the equation 
sls = dsolve(eq1) # solvde ODE 

# print solutions: 
print("For ode") 
display(eq1) 
print("the solutions are:") 
for s in sls: 
    display(s) 

# plot solutions: 
x = np.linspace(0, 2, 100) 
fg, axx = plt.subplots(2, 1) 
axx[0].set_title("Real part of solution of $\\frac{d}{dt}f(t)= (f(t))^4$") 
axx[1].set_title("Imag. part of solution of $\\frac{d}{dt}f(t)= (f(t))^4$") 
fg.suptitle("$C_1=0.1$") 
for i, s in enumerate(sls, start=1): 
    fn1 = s.rhs.subs("C1", .1) # C_1 -> 1 
    fn2 = sy.lambdify(t, fn1, modules="numpy") # make numpy function 
    y = fn2(x+0j) # needs to be called with complex number 
    axx[0].plot(x, np.real(y), label="Sol. %d" % i) 
    axx[1].plot(x, np.imag(y), label="Sol. %d" % i) 
for ax in axx: 
    ax.legend(loc="best") 
    ax.grid(True) 
axx[0].set_ylabel("Re$\\{f(t)\\}$") 
axx[1].set_ylabel("Im$\\{f(t)\\}$") 
axx[-1].set_xlabel("$t$") 
fg.canvas.draw() 
plt.show() 

В оболочке IPython вы должны увидеть следующее: