【问题标题】:Solving nonlinear differential first order equations using Python使用 Python 求解非线性微分一阶方程
【发布时间】:2015-08-21 15:48:36
【问题描述】:

我想使用 Python 求解非线性一阶微分方程。

例如,

df/dt = f**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'

【问题讨论】:

  • 试试t=np.linspace(0,0.00033,1000)
  • 这很奇怪,因为当我将scipy.integrate.odeint 更改为odeint 时,上面的代码对我来说工作正常(因为您实际上没有导入 scipy 命名空间,所以应该只按其名称调用 odeint 函数)
  • @HYRY :我试过t=np.linspace(0,0.00033,1000),它适用于这个t,我怎样才能让它适用于t=np.linspace(0,2,100)? @Azrathud:它符合您的建议......但并非一直如此。实际上,如果你尝试我的代码 10 次,它不会工作 10 次:会有AttributeError。你有其他建议让它发挥作用吗?
  • @Jack,请查看符号解决方案:wolframalpha.com/input/?i=df%2Fdt%20%3D%20f**4
  • @HYRY :我将您的链接与df/dt = f^4, f(0)=10 一起使用,WolframAlpha 给出了f(x) = 10 / (1 - 3000x)^(1 / 3) 作为解决方案。然后我在 Geogebra 上绘制了解决方案,它与我运行代码后绘制的结果不一样。

标签: python math numpy matplotlib scipy


【解决方案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 shell中,您应该看到以下内容:

【讨论】:

  • 谢谢你,非常有趣,我已经尝试了一个复杂的非线性方程式的解决方案,但它永远保持嘎吱作响,有没有办法限制计算时间,或者我错过了什么,等等:eq1 = sy.eq(t *(3 * f(t).diff(t)-f(t).diff(t)** 3),f(t)*(1-3 * f(t).diff (t)** 2))和dsolve:sls = dsolve(eq1,ics = {f(0):-2,f(t).diff(t).subs(t,0):0})
  • @ fabianotarlao不确定是否为您的等式存在封闭式解决方案 - 乍一看它似乎并不是那么... span>
猜你喜欢
  • 1970-01-01
  • 2021-01-29
  • 2021-06-07
  • 1970-01-01
  • 2021-06-14
  • 2020-08-18
  • 1970-01-01
  • 1970-01-01
  • 1970-01-01
相关资源
最近更新 更多