Python中数值类型:Python中的数值 ODE求解

如何在 Python 中数值求解 ODE?

考虑

equation to solve

\ddot{u}(\phi) = -u + \sqrt{u}

具有以下条件

u(0) = 1.49907

\dot{u}(0) = 0

有了约束

0 <= \phi <= 7\pi.

最后,我想产生一个参数图,其中 x 和 y 坐标作为 u 的函数生成。

问题是,我需要运行 odeint 两次,因为这是一个二阶微分方程。我试图让它在第一次之后再次运行,但它返回一个雅可比错误。必须有一种方法一次运行两次。

这里是错误:

odepack.error:函数及其 Jacobian 必须是可调用函数

的代码生成。有问题的行是 sol = odeint。

import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt
from numpy import linspace
def f(u, t):
    return -u + np.sqrt(u)
times = linspace(0.0001, 7 * np.pi, 1000)
y0 = 1.49907
yprime0 = 0
yvals = odeint(f, yprime0, times)
sol = odeint(yvals, y0, times)
x = 1 / sol * np.cos(times)
y = 1 / sol * np.sin(times)
plot(x,y)
plt.show()
Edit

我正在尝试构建第 9 页的情节

Classical Mechanics Taylor

这是 Mathematica 的情节

mathematica plot

In[27]:= sol = 
 NDSolve[{y''[t] == -y[t] + Sqrt[y[t]], y[0] == 1/.66707928, 
   y'[0] == 0}, y, {t, 0, 10*\[Pi]}];
In[28]:= ysol = y[t] /. sol[[1]];
In[30]:= ParametricPlot[{1/ysol*Cos[t], 1/ysol*Sin[t]}, {t, 0, 
  7 \[Pi]}, PlotRange -> {{-2, 2}, {-2.5, 2.5}}]
27
import scipy.integrate as integrate
import matplotlib.pyplot as plt
import numpy as np
pi = np.pi
sqrt = np.sqrt
cos = np.cos
sin = np.sin
def deriv_z(z, phi):
    u, udot = z
    return [udot, -u + sqrt(u)]
phi = np.linspace(0, 7.0*pi, 2000)
zinit = [1.49907, 0]
z = integrate.odeint(deriv_z, zinit, phi)
u, udot = z.T
# plt.plot(phi, u)
fig, ax = plt.subplots()
ax.plot(1/u*cos(phi), 1/u*sin(phi))
ax.set_aspect('equal')
plt.grid(True)
plt.show()

enter image description here

4

您的其他question的代码非常接近您想要的。需要进行两个更改:

你正在解决一个不同的 ODE(因为你改变了函数deriv中的两个符号)

所需绘图的y组件来自解值,而不是来自解的一阶导数的值,因此您需要将u[:,0](函数值)替换为u[:, 1](导数)。

这是最终结果:

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint
def deriv(u, t):
    return np.array([u[1], -u[0] + np.sqrt(u[0])])
time = np.arange(0.01, 7 * np.pi, 0.0001)
uinit = np.array([1.49907, 0])
u = odeint(deriv, uinit, time)
x = 1 / u[:, 0] * np.cos(time)
y = 1 / u[:, 0] * np.sin(time)
plt.plot(x, y)
plt.show()

但是,我建议您使用 unutbu 的答案中的代码,因为它是自记录(u, udot = z),并使用np.linspace而不是np.arange

x = 1 / u * np.cos(phi)
y = 1 / u * np.sin(phi)
plt.plot(x, y)
plt.show()
3

您可以使用 scipy.integrate.ode。要解决 dy / dt = f(t,y),初始条件 y(t0)= y0,在时间 = t1 与四阶 Runge-Kutta 你可以做这样的事情:

from scipy.integrate import ode
solver = ode(f).set_integrator('dopri5')
solver.set_initial_value(y0, t0)
dt = 0.1
while t < t1:
    y = solver.integrate(t+dt)
    t += dt

编辑:你必须让你的导数一阶使用数值积分。这可以通过设置例如 z1 = u 和 z2 = du / dt 来实现,之后你有 dz1 / dt = z2 和 dz2 / dt = d ^ 2u / dt ^ 2。将这些代入你的原始方程,并简单地迭代向量 dZ / dt,这是一阶。

编辑 2:这是整个事情的示例代码:

import numpy as np
import matplotlib.pyplot as plt
from numpy import sqrt, pi, sin, cos
from scipy.integrate import ode
# use z = [z1, z2] = [u, u']
# and then f = z' = [u', u''] = [z2, -z1+sqrt(z1)]
def f(phi, z):
    return [z[1], -z[0]+sqrt(z[0])]
# initialize the 4th order Runge-Kutta solver
solver = ode(f).set_integrator('dopri5')
# initial value
z0 = [1.49907, 0.]
solver.set_initial_value(z0)
values = 1000
phi = np.linspace(0.0001, 7.*pi, values)
u = np.zeros(values)
for ii in range(values):
    u[ii] = solver.integrate(phi[ii])[0] #z[0]=u
x = 1. / u * cos(phi)
y = 1. / u * sin(phi)
plt.figure()
plt.plot(x,y)
plt.grid()
plt.show()
2

scipy.integrate()进行 ODE 集成。这就是您要找的吗?

本站系公益性非盈利分享网址,本文来自用户投稿,不代表边看边学立场,如若转载,请注明出处

(729)
注册表cmd:私有npm注册表 回退到全局注册表
上一篇
门机变频器怎么调试:环形振荡器改变频率(ring oscillator)
下一篇

相关推荐

  • python精度:如何利用Python来提高精度

    Python精度是指Python程序在数值计算时所能提供的最大精度。它取决于Python程序使用的数字类型,以及Python解释器的精度。…

    2023-04-01 04:27:43
    0 22 89
  • python界面开发构建一个简单、可靠的用户界面

    Python界面开发是指使用Python语言来创建图形用户界面(GUI)的过程。它可以帮助你创建可视化的应用程序,使用户能够与你的程序交互。…

    2023-03-18 10:39:36
    0 88 29
  • python 软件测试深入了解如何确保质量

    Python软件测试是一种使用Python编写的自动化测试,它可以帮助开发人员进行软件测试。Python软件测试主要涉及到三个方面:单元测试、集成测试和系统测试。…

    2023-08-14 11:34:51
    0 35 57
  • pythonlist切片:利用Python列表切片获取最大价值

    示例示例Python列表切片是指从列表中提取特定范围的元素,而不需要遍历整个列表。 列表切片使用方括号[]和冒号:来表示,其形式如下:list[start:end:step]…

    2023-07-31 01:45:29
    0 59 47
  • python imread函数一步搞定!

    Python imread函数是用于从图像文件读取图像数据的函数,它是由Scipy库提供的。它可以从多种格式的图像文件中读取图像数据,包括JPEG、PNG、TIFF、GIF等。…

    2023-05-26 15:06:44
    0 56 47
  • python 创建列表:使用Python创建一个强大的列表

    示例示例Python创建列表的方法有多种,下面介绍其中几种常用的方法。使用 [] 创建空列表…

    2023-03-03 12:56:32
    0 98 47
  • pythonweb服务器:如何使用Python搭建Web服务器

    Python Web服务器是一种使用Python语言编写的Web服务器,它可以接受HTTP请求并返回相应的响应。Python Web服务器可以处理动态内容,比如数据库查询,文件上传,CGI脚本等。…

    2023-01-24 10:13:27
    0 50 98
  • python窗口代码从入门到精通

    Python窗口代码是使用Python语言创建GUI(图形用户界面)应用程序的一种方式。它使用Python的tkinter模块,可以快速创建简单的窗口,并使用其中的控件进行交互。下面是一个简单的窗口示例:…

    2023-07-04 11:48:45
    0 32 19

发表评论

登录 后才能评论

评论列表(10条)