求大神看看哪错了,怎么改
oqvpo 2019-04-22 07:45:27 import numpy as np
import matplotlib.pyplot as plt
from math import *
time=eval(input('time=?'))
dt=0.01
n=int(round(time/dt))
t=np.zeros(n,float)
E=np.zeros(n,float)
w=np.zeros(n,float)
theta=np.zeros(n,float)
x=np.zeros(n,float)
T=np.zeros(n,float)
V=np.zeros(n,float)
P=np.zeros(n,float)
Ah=np.zeros(n,float)
Al=np.zeros(n,float)
p=1.2
Vc=40*10e-6
cV=717.0
Th=600.0
Tl=300.0
u=100000.0
A0=40*10e-4
Ap=1.9*10e-4
A1_max=4.9*10e-2
R=1.25*10e-2
I=4*10e-4
b=0.7*10e-4
P0=10**5
m=p*(Vc+Ap*R)
P[0]=P0
E[0]=cV*m*300
w[0]=0.0
theta[0]=0.0
t[0]=0.0
x[0]=R
V[0]=Vc+Ap*x[0]
T[0]=560
for i in range(n):
Al[i]=A0*(1-cos(theta[i]))/2+A1_max*x[i]
Ah[i]=A0*(1+cos(theta[i]))/2
E[i+1]=E[i]+Ah[i]*u*(Th-T[i])-(Al[i]*u*(T[i]-Tl)-(P[i]-P0)*Ap*R*w[i]*cos(theta[i]))*dt
w[i+1]=w[i]+((P[i]-P0)*Ap*R*cos(theta[i])-b*w[i])/I*dt
theta[i+1]=theta[i]+w[i]*dt
print(E[i])
x[i+1]=R*(1+sin(theta[i+1]))
T[i+1]=E[i+1]/(cV*m)
V[i+1]=Vc+Ap*x[i+1]
P[i+1]=2/5*E[i+1]/V[i+1]
t[i+1]=t[i]+dt
plt.plot(t,x)
plt.xlabel('t')
plt.ylabel('x')
plt.show()
输出:time=?1
164.5515
144239.5515
-1065418312.7497553
7802303378666.067
3.590169491968912e+18
-1.305296690914636e+28
2.754751549541168e+43
8.04321363123899e+67
-2.6029319550428185e+107
3.7708738042562866e+171
Warning (from warnings module):
File "C:\Users\Administrator\Desktop\Project striling engine.py", line 42
E[i+1]=E[i]+Ah[i]*u*(Th-T[i])-(Al[i]*u*(T[i]-Tl)-(P[i]-P0)*Ap*R*w[i]*cos(theta[i]))*dt
RuntimeWarning: overflow encountered in double_scalars
1.1942494159483741e+275
Warning (from warnings module):
File "C:\Users\Administrator\Desktop\Project striling engine.py", line 42
E[i+1]=E[i]+Ah[i]*u*(Th-T[i])-(Al[i]*u*(T[i]-Tl)-(P[i]-P0)*Ap*R*w[i]*cos(theta[i]))*dt
RuntimeWarning: invalid value encountered in double_scalars
-inf
nan
Traceback (most recent call last):
File "C:\Users\Administrator\Desktop\Project striling engine.py", line 46, in <module>
x[i+1]=R*(1+sin(theta[i+1]))
ValueError: math domain error