2015-01-07 38 views
5

私は与えられた初期条件で、この微分方程式を解決したい:Pythonの組み込み関数odeintを使って微分方程式を解く方法は?

(3x-1)y''-(3x+2)y'+(6x-8)y=0, y(0)=2, y'(0)=3 

ANSがここ

y=2*exp(2*x)-x*exp(-x)

でなければなりませんが、私のコードです:

def g(y,x): 
    y0 = y[0] 
    y1 = y[1] 
    y2 = (6*x-8)*y0/(3*x-1)+(3*x+2)*y1/(3*x-1) 
    return [y1,y2] 

init = [2.0, 3.0] 
x=np.linspace(-2,2,100) 
sol=spi.odeint(g,init,x) 
plt.plot(x,sol[:,0]) 
plt.show() 

しかし、私は取得します答えとは異なります。 どうしたのですか?

答えて

12

ここにはいくつかの問題があります。まず、あなたの方程式は明らかに明らかです。

(3x-1)y '' - (3x + 2)y ' - (6x-8)y = 0; y(0)= 2、y '(0)= 3

(yの項の記号に注意)。この式では、解析解とy2の定義が正しいです。 '' 'Y、及びそれらの誘導体を返すとy y[0](Y)、y[1](Y)':@Warren Weckesserが言うよう

第二に、あなたはgからyとして2つのパラメータを渡す必要があります。

第3に、初期条件はx = 0ですが、積分するxグリッドは-2で始まります。 odeintためのドキュメントから、このパラメータを、それらのコールサインの説明でt

odeint(func, y0, t, args=(),...)

T:yについて解くための時間ポイントの配列 配列。最初の 値ポイントは、このシーケンスの最初の要素である必要があります。

したがって、0から開始して統合するか、-2から始まる初期条件を提供する必要があります。

最後に、統合の範囲はx = 1/3で特異点をカバーします。 odeintはここでは時間がかかるかもしれませんが(明らかにそうではありません)

ここで働くようで一つのアプローチです:二階微分方程式については

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

def g(y, x): 
    y0 = y[0] 
    y1 = y[1] 
    y2 = ((3*x+2)*y1 + (6*x-8)*y0)/(3*x-1) 
    return y1, y2 

# Initial conditions on y, y' at x=0 
init = 2.0, 3.0 
# First integrate from 0 to 2 
x = np.linspace(0,2,100) 
sol=odeint(g, init, x) 
# Then integrate from 0 to -2 
plt.plot(x, sol[:,0], color='b') 
x = np.linspace(0,-2,100) 
sol=odeint(g, init, x) 
plt.plot(x, sol[:,0], color='b') 

# The analytical answer in red dots 
exact_x = np.linspace(-2,2,10) 
exact_y = 2*np.exp(2*exact_x)-exact_x*np.exp(-exact_x) 
plt.plot(exact_x,exact_y, 'o', color='r', label='exact') 
plt.legend() 

plt.show() 

enter image description here

+1

は、 'init'は、長さ2を持っていない3はずです(と' G'は長さ2を返す必要がありますアレイ)。 –

+0

あなたはそうです:私は混乱しています。私はそれを修正するために編集しました。 – xnx

関連する問題