我正在求解一个耦合微分方程组，其中一个"；常数"；在微分方程中，实际上是一个周期性变化的值：周期的前半部分为1，其余部分为0，周期为2pi。
import matplotlib.pyplot as plt
import numpy as np
import math
l=1
lt=0
k=0.05
eps= 0.11
w= (1 - eps)
w_0= (1 + eps)
u=(w_0/4*w)*(w**2+(k/2)**2)
def RungeKutta(t, jx, jy, jz, x, p, dt, djxdt, djydt, djzdt, dxdt, dpdt):
a1 = dt * djxdt(t, jx, jy, jz, x, p)
b1 = dt * djydt(t, jx, jy, jz, x, p)
c1 = dt * djzdt(t, jx, jy, jz, x, p)
d1 = dt * dxdt(t, jx, jy, jz, x, p)
e1 = dt * dpdt(t, jx, jy, jz, x, p)
a2 = dt * djxdt(t + 0.5 * dt, jx + 0.5 * a1, jy + 0.5 * b1, jz + 0.5 * c1, x + 0.5 * d1, p + 0.5 * e1)
b2 = dt * djydt(t + 0.5 * dt, jx + 0.5 * a1, jy + 0.5 * b1, jz + 0.5 * c1, x + 0.5 * d1, p + 0.5 * e1)
c2 = dt * djzdt(t0 + 0.5 * dt, jx + 0.5 * a1, jy + 0.5 * b1, jz + 0.5 * c1, x + 0.5 * d1, p + 0.5 * e1)
d2 = dt * dxdt(t0 + 0.5 * dt, jx + 0.5 * a1, jy + 0.5 * b1, jz + 0.5 * c1, x + 0.5 * d1, p + 0.5 * e1)
e2 = dt * dpdt(t0 + 0.5 * dt, jx + 0.5 * a1, jy + 0.5 * b1, jz + 0.5 * c1, x + 0.5 * d1, p + 0.5 * e1)
a3 = dt * djxdt(t + 0.5 * dt, jx + 0.5 * a2, jy + 0.5 * b2, jz + 0.5 * c2, x + 0.5 * d2, p + 0.5 * e2)
b3 = dt * djydt(t + 0.5 * dt, jx + 0.5 * a2, jy + 0.5 * b2, jz + 0.5 * c2, x + 0.5 * d2, p + 0.5 * e2)
c3 = dt * djzdt(t + 0.5 * dt, jx + 0.5 * a2, jy + 0.5 * b2, jz + 0.5 * c2, x + 0.5 * d2, p + 0.5 * e2)
d3 = dt * dxdt(t + 0.5 * dt, jx + 0.5 * a2, jy + 0.5 * b2, jz + 0.5 * c2, x + 0.5 * d2, p + 0.5 * e2)
e3 = dt * dpdt(t + 0.5 * dt, jx + 0.5 * a2, jy + 0.5 * b2, jz + 0.5 * c2, x + 0.5 * d2, p + 0.5 * e2)
a4 = dt * djxdt(t + dt, jx + a3, jy + b3, jz + c3, x + d3, p + e3)
b4 = dt * djydt(t + dt, jx + a3, jy + b3, jz + c3, x + d3, p + e3)
c4 = dt * djzdt(t + dt, jx + a3, jy + b3, jz + c3, x + d3, p + e3)
d4 = dt * dxdt(t + dt, jx + a3, jy + b3, jz + c3, x + d3, p + e3)
e4 = dt * dpdt(t + dt, jx + a3, jy + b3, jz + c3, x + d3, p + e3)
t = t+dt
jx = jx+(1/6)*(a1+2*a2+2*a3+a4)
jy = jy+(1/6)*(b1+2*b2+2*b3+b4)
jz = jz+(1/6)*(c1+2*c2+2*c3+c4)
x = x+(1/6)*(d1+2*d2+2*d3+d4)
p = p+(1/6)*(e1+2*e2+2*e3+e4)
return t, jx, jy, jz, x, p
def djxdt(t, jx, jy, jz, x, p):
return -w_0*jy
def djydt(t, jx, jy, jz, x, p):
return w_0*jx-2*lambda_t(t)*(math.sqrt(2*w))*x*jz
def djzdt(t, jx, jy, jz, x, p):
return 2*lambda_t(t)*(math.sqrt(2*w))*x*jy
def dxdt(t, jx, jy, jz, x, p):
return p-(k/2)*x
def dpdt(t, jx, jy, jz, x, p):
return -(w**2)*x-(k/2)*p-2*lambda_t(t)*(math.sqrt(2*w))*jx
def lambda_t(t):
if math.fmod(t, 2*np.pi) > np.pi:
return 1
else:
return 0
t0=0
jx0= 0.5*math.sqrt(1-u**2)
jy0= 0
jz0= -0.5*u
x0= -math.sqrt(2*w*(1-u**2))/(w**2-(k**2)/4)
p0= -(k/2)*math.sqrt(2*w*(1-u**2))/(w**2-(k**2)/4)
dt=0.0001*2*np.pi
t_end=50
t_list = [t0]
jx_list = [jx0]
jy_list = [jy0]
jz_list = [jz0]
x_list = [x0]
p_list = [p0]
t_list2 = [t0]
jx_list2 = [jx0]
t = t0
jx=jx0
jy=jy0
jz=jx0
x=x0
p=p0
while t <= t_end:
t, jx, jy, jz, x, p = RungeKutta(t, jx, jy, jz, x, p, dt, djxdt, djydt, djzdt, dxdt, dpdt)
if (t%2*math.pi<0.000001 or t%2*math.pi>0.000001):
t_list2.append(t)
jx_list2.append(jx)
t_list.append(t)
jx_list.append(jx)
jy_list.append(jy)
jz_list.append(jz)
x_list.append(x)
p_list.append(p)
plt.plot(t_list, jx_list, label="Jx")
plt.plot(t_list, jy_list, label="Jy")
plt.plot(t_list, jz_list, label="Jz")
plt.legend(fontsize=15)
plt.show()
plt.scatter(t_list2, jx_list2)
plt.show()
我不知道为什么总是出现未缩进的错误与任何外部缩进级别都不匹配。
此外，我只需要在t是2pi的倍数时绘制点和值。我试过的代码
if (t%2*math.pi==0):
t_list2.append(t)
jx_list2.append(jx)
但这个数字在t=0时只有一个点如何解决此问题？我以前的设置会影响这一部分吗？
谢谢你的帮助！
从纯python的角度来看，这个函数中有几个错误
def l(f):
if int(2*t)%2np.pi == 0:
return 1
else:
return -1
- 函数中未使用参数
f，而是显示未声明的t 2np.pi应该立即出现语法错误int(2*t)%2*np.pi是(int(2*t)%2)*np.pi，这可能不是期望的解释。使用int((2*f)%(2*np.pi))。您也可以取消2- 较低的值为
-1，而在描述中则为0
model是一个由3个状态组件组成的系统，您的初始状态有5个组件。这是不兼容的
- 应将重复表达式分配给局部变量，以避免重复计算
请报告详细的错误描述，而不仅仅是您对观察到的错误的总结。