2017-03-06 13 views
1

私は式dy/dx = x + y/5と初期値y(0) = -3を持っています。numpyで一次微分方程式を表現する

pyplotを使ってこの関数の正確なグラフをプロットする方法を知りたいと思います。

私はx軸として使用したいx = np.linspace(0, interval, steps+1)も持っています。だから私はy軸の値だけを探しています。

ありがとうございます。

+0

。 –

+0

これが役立つかどうかを確認してください:http://stackoverflow.com/questions/16001341/having-trouble-while-using-scipy-integrate-odeint-with-python –

答えて

2

d.e.を解決する必要があるとすれば、あなたはsympyを使って代数的にこれをやってほしいかもしれません。 (そうでないかもしれません)

モジュールをインポートし、関数と従属変数を定義します。

>>> from sympy import * 
>>> f = Function('f') 
>>> var('x') 
x 

ソルバを起動します。 d.e。等号の左に転記する必要があり、yをその関数の指定子に置き換える必要があります。

>>> dsolve(Derivative(f(x),x)-x-f(x)/5) 
Eq(f(x), (C1 + 5*(-x - 5)*exp(-x/5))*exp(x/5)) 

期待どおり、解は任意の定数で与えられます。初期値を使って解決する必要があります。これをsympy変数と定義します。

​​

ここで、この任意の定数を、私たちが解くことができる方程式の左側として表す式を作成します。 f(0)を初期状態の値に置き換えます。次に、この条件でxの値を代入して、C1の式を取得します。

>>> expr = -3 - ((C1 + 5*(-x - 5)*exp(-x/5))*exp(x/5)) 
>>> expr.subs(x,0) 
-C1 + 22 

つまり、C1 = 22です。最後に、この値を使用して微分方程式の特定の解を得ることができます。

>>> ((C1 + 5*(-x - 5)*exp(-x/5))*exp(x/5)).subs(C1,22) 
((-5*x - 25)*exp(-x/5) + 22)*exp(x/5) 

私はぼんやりと私は、この関数は、初期条件を満たしていることを確認ひどい間違いを、これまで恐ろしいだから。

>>> (((-5*x - 25)*exp(-x/5) + 22)*exp(x/5)).subs(x,0) 
-3 

(通常物事が間違っている、私はそれらをチェックするのを忘れた場合にのみ。このような人生がある。)

そして、私はあまりにもsympyでこれをプロットすることができます。

>>> plot(((-5*x - 25)*exp(-x/5) + 22)*exp(x/5),(x,-1,5)) 
<sympy.plotting.plot.Plot object at 0x0000000008C2F780> 

solution to d.e.

3

ただ、完全を期すため、式のこの種のは簡単scipy.integrate.odeintを使用して、数値的に統合することができます。あなたが最初にあなたの微分方程式を解くために必要

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

# function dy/dx = x + y/5. 
func = lambda y,x : x + y/5. 
# Initial condition 
y0 = -3 # at x=0 
# values at which to compute the solution (needs to start at x=0) 
x = np.linspace(0, 4, 101) 
# solution 
y = odeint(func, y0, x) 
# plot the solution, note that y is a column vector 
plt.plot(x, y[:,0]) 
plt.xlabel('x') 
plt.ylabel('y') 
plt.show() 

enter image description here

関連する問題