funcion - odeint python install

Resolución numérica de ODE en Python

Question

¿Cómo puedo resolver numéricamente una ODE en Python?

Considerar

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

con las siguientes condiciones

u(0) = 1.49907

y

/dot{u}(0) = 0

con la restricción

0 <= /phi <= 7/pi.

Entonces, finalmente, quiero producir una gráfica paramétrica donde las coordenadas xey se generan como una función de u.

El problema es que necesito ejecutar odeint dos veces ya que esta es una ecuación diferencial de segundo orden. Intenté que se ejecutara nuevamente después de la primera vez, pero vuelve con un error jacobiano. Debe haber una forma de ejecutarlo dos veces a la vez.

Aquí está el error:

odepack.error: la función y su Jacobian deben ser funciones invocables

que genera el código a continuación. La línea en cuestión es el 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()

Editar

Estoy tratando de construir la trama en la página 9

Mecánica clásica Taylor

Aquí está la trama con Mathematica

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}}]

Answer 1

scipy.integrate () integra ODE. ¿Es esto lo que estás buscando?

Answer 2

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()

Answer 3

Puedes usar scipy.integrate.ode. Para resolver dy / dt = f (t, y), con la condición inicial y (t0) = y0, at time = t1 con Runge-Kutta de 4to orden, podrías hacer algo como esto:

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

Editar: Tienes que obtener tu derivada al primer orden para usar la integración numérica. Esto lo puede lograr estableciendo, por ejemplo, z1 = u y z2 = du / dt, después de lo cual tiene dz1 / dt = z2 y dz2 / dt = d ^ 2u / dt ^ 2. Sustituya estos en su ecuación original, y simplemente itere sobre el vector dZ / dt, que es de primer orden.

Editar 2: Aquí hay un código de ejemplo para todo:

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()

Answer 4

El código de tu otra pregunta está muy cerca de lo que quieres. Se necesitan dos cambios:

• Estabas resolviendo una EDO diferente (porque cambiaste dos signos dentro de la función deriv )
• El componente y de su trazado deseado proviene de los valores de solución, no de los valores de la primera derivada de la solución, por lo que debe reemplazar u[:,0] (valores de función) por u[:, 1] (derivadas) .

Este es el resultado final:

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()

Sin embargo, le sugiero que use el código de la respuesta de unutbu porque se auto documenta ( u, udot = z ) y usa np.linspace lugar de np.arange . Luego, ejecuta esto para obtener tu figura deseada:

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