用python解決二階ODE

我在寫一些代碼來求解二階微分方程，但是它給出了一個完全錯誤的結果。我認爲問題在於以正確的方式表達歐拉方法。我也嘗試了另一個二階ODE，但我也未能近似y（x）。你能指出哪裏可能出錯嗎？請參閱圖形和代碼：用python解決二階ODE

求解ODE：

y"(x)=-1001y'(x)-1000y(x) 

y(0)=1, y'(0)=0

解析解：

y(x)=(1000exp(-x)-exp(-1000x))/999

重寫爲2個常微分方程的系統;

y'(x)=v(x) 

v'(x)=-1001v(x)-1000y(x)=f(x,y,v) 

y(0)=1,v(0)=0

在Jupyter筆記本的代碼（蟒蛇）：

import numpy as np 
import matplotlib.pyplot as plt 
%matplotlib inline 

a=-0.0020 
b=0.0 
h=0.0005 #step size 
x=np.arange(a,b,h) 
n=int((b-a)/h) 

def y_exact(x): 
    return (1000*np.exp(-x)-np.exp(-1000*x))/999 
def f(y,v): 
    return -1001*v-1000*y 

y_eu=np.zeros(n,float) 
y_ei=np.zeros(n,float) 

y_eu[0]=1.0 
y_ei[0]=1.0 

v_eu[0]=0.0 
v_ei[0]=0.0 

for j in range(0,n-1): 
    #Euler's method (explicit) 
    y_eu[j+1]=y_eu[j]+v_eu[j]*h 
    v_eu[j+1]=v_eu[j]+f(y_eu[j],v_eu[j])*h 

    #Implicit Euler's method 
    v_ei[j+1]=(v_ei[j]-1000*y_ei[j]*h)/(1+1001*h+1000*h**2) 
    y_ei[j+1]=y_ei[j]+v_ei[j+1]*h 
plt.plot(x,y_exact(x),label="Exact solution") 
plt.plot(x,y_eu,label="Euler's method") 
plt.plot(x,y_ei,label="Implicit Euler's method") 
plt.legend() 
plt.show()

謝謝！

來源

2017-04-25 user7918565

爲什麼試圖再次實現它們？你可以使用[sympy]（http://docs.sympy.org/latest/modules/solvers/ode.html#ode-docs） –

和你的'y_eu [0]'等不正確。 'y_eu [0]'是-0.0020的值，而不是0 –

或者換句話說，爲了反映使用的初始值和積分方向，你必須使用'plot（x，y_exact（xa），label = 「精確的解決方案」），當將「h」減少10倍時，它給出了很好的視覺融合。 – LutzL

感謝評論者Maarten Fabre的回答。代碼無法正常工作，因爲初始條件，即歐拉方法的y0和v0被錯誤地採用：

def y_exact(x): 
    return (1000*np.exp(-x)-np.exp(-1000*x))/999 
def dydx_exact(x): 
    return -1000*(np.exp(-x)-np.exp(-1000*x))/999 

y_eu[0]=y_exact(a) 
v_eu[0]=dydx_exact(a)

來源

2017-04-25 16:53:13 user7918565

用python解決二階ODE

回答

相關問題