2016-10-16 45 views
1

我想從這個網站了解它:http://nbviewer.jupyter.org/github/numerical-mooc/numerical-mooc/blob/master/lessons/01_phugoid/01_03_PhugoidFullModel.ipynb試圖解決數值DIFF式使用歐拉法,無效值錯誤

我試圖把它用盡可能少的幫助,可能的代碼,但我一直收到此錯誤:

C:\用戶\ 「我的真名」 \ Anaconda2 \ LIB \站點包\ ipykernel__main __潘岳:29:RuntimeWarning:在double_scalars

隨着我的陰謀沒有數據點遇到無效值。所以我直接從網站上直接粘貼了所有的代碼,我仍然遇到錯誤!我放棄了,有人可以幫助一個蟒蛇新手嗎?

import numpy as np 
from matplotlib import pyplot 
from math import sin, cos, log, ceil 
%matplotlib inline 
from matplotlib import rcParams 
rcParams['font.family'] = 'serif' 
rcParams['font.size'] = 16 

# model parameters: 
g = 9.8  # gravity in m s^{-2} 
v_t = 30.0 # trim velocity in m s^{-1} 
C_D = 1/40 # drag coefficient --- or D/L if C_L=1 
C_L = 1 # for convenience, use C_L = 1 

### set initial conditions ### 
v0 = v_t  # start at the trim velocity (or add a delta) 
theta0 = 0 # initial angle of trajectory 
x0 = 0  # horizotal position is arbitrary 
y0 = 1000 # initial altitude 



def f(u): 

    v = u[0] 
    theta = u[1] 
    x = u[2] 
    y = u[3] 
    return np.array([-g*sin(theta) - C_D/C_L*g/v_t**2*v**2, -g*cos(theta)/v + g/v_t**2*v, v*cos(theta), v*sin(theta)]) 

def euler_step(u, f, dt): 
    u + dt * f(u) 

T = 100       # final time 
dt = 0.1       # time increment 
N = int(T/dt) + 1     # number of time-steps 
t = np.linspace(0, T, N)  # time discretization 

# initialize the array containing the solution for each time-step 
u = np.empty((N, 4)) 
u[0] = np.array([v0, theta0, x0, y0])# fill 1st element with initial values 

# time loop - Euler method 
for n in range(N-1): 
    u[n+1] = euler_step(u[n], f, dt) 





x = u[:,2] 
y = u[:,3] 


pyplot.figure(figsize=(8,6)) 
pyplot.grid(True) 
pyplot.xlabel(r'x', fontsize=18) 
pyplot.ylabel(r'y', fontsize=18) 
pyplot.title('Glider trajectory, flight time = %.2f' % T, fontsize=18) 
pyplot.plot(x,y, 'k-', lw=2); 

回答

0

的解決方案是非常簡單的。你忘記了euler_step中的return語句。 變化

def euler_step(u, f, dt): 
    u + dt * f(u) 

def euler_step(u, f, dt): 
    return u + dt * f(u) 

,也將努力

+0

謝謝你,但現在我得到一個平坦的直線爲我的y值!你可以運行代碼,看看你是否得到相同的? –

+0

我得到了一些向下彎曲的軌跡,所以它似乎工作。你可以直接看看y值並檢查它是一個繪圖還是計算問題? – Jannick

+0

很奇怪,我打印了我的y值,它們都是1000.所有這些都是恆定的。 –