2016-10-27 6 views
0

を間違って何をやっている私は、Pythonに完全に新しいですし、次の頌歌を統合しよう:私はこのDopri5実装で

$ \ドット{X} = -2x-y^2 $

$ \ドット{y} = -yx^2

この結果、すべて0の配列になります。 私は間違っていますか?それは主にコピーされたコードであり、別のものと結合されたコードではうまく機能しませんでした。

import numpy as np 
import matplotlib.pyplot as plt 
from scipy.integrate import ode 


def fun(t, z): 
    """ 
    Right hand side of the differential equations 
     dx/dt = -omega * y 
     dy/dt = omega * x 
    """ 
    x, y = z 
    f = [-2*x-y**2, -y-x**2] 
    return f 

# Create an `ode` instance to solve the system of differential 
# equations defined by `fun`, and set the solver method to 'dop853'. 
solver = ode(fun) 
solver.set_integrator('dopri5') 

# Set the initial value z(0) = z0. 
t0 = 0.0 
z0 = [0, 0] 
solver.set_initial_value(z0, t0) 

# Create the array `t` of time values at which to compute 
# the solution, and create an array to hold the solution. 
# Put the initial value in the solution array. 
t1 = 2.5 
N = 75 
t = np.linspace(t0, t1, N) 
sol = np.empty((N, 2)) 
sol[0] = z0 

# Repeatedly call the `integrate` method to advance the 
# solution to time t[k], and save the solution in sol[k]. 
k = 1 
while solver.successful() and solver.t < t1: 
    solver.integrate(t[k]) 
    sol[k] = solver.y 
    k += 1 

# Plot the solution... 
plt.plot(t, sol[:,0], label='x') 
plt.plot(t, sol[:,1], label='y') 
plt.xlabel('t') 
plt.grid(True) 
plt.legend() 
plt.show() 

答えて

1

あなたの初期状態(z0)が[0,0]です。この初期状態の時間導関数(fun)も[0,0]です。したがって、この初期条件では、常に[0,0]が正しい解です。

初期条件を他の値に変更すると、より興味深い結果が得られます。

+0

ああ、それは固定小数点であり、私はコード内でsomethignをやっていると思った –

関連する問題