我想在python代码中使用Dirac delta作为时间函数来解决四个耦合问题
微分方程。在python代码中,我使用solve_ivp来求解耦合方程。
我没弄对。我无法定义狄拉克δ函数。如果需要,请帮助我
任何人都可以。在python代码中,我使用solve_ivp来求解耦合方程。
我没弄对。我无法定义狄拉克δ函数。如果需要,请帮助我
任何人都可以
import numpy as np
import matplotlib.pyplot as plt
import scipy.interpolate
from pylab import *
from qutip import *
import scipy.special as sp
import scipy.linalg as la
from scipy.integrate import solve_ivp
import math
import cmath
from sympy import DiracDelta, diff, pi
from scipy import signal
omg1 = 1
omg2 = 1.1
Omg=1
Omegar=10
beta=0.1
Mu=0.5
epsilon=0.01
V=1
K=0.5
g=1
m=1
hcut =1
A = 0.001
p = [ omg1, omg2, Omg,Omegar, beta, Mu, K, g,m]
#------------------------------------------------------------------------------
#####----------------Initial conditions, packed in w0--------------------------
##### IMPORTANT NOTE:
##### Please feed initial values with a complex part even if it's zero
#y1 & y2 are first derivatives of x1 and x2
x1 = 1
x2 = 0
#------------------------------------------------------------------------------
z1 = 0+0j
z2 = 1+0j
#A = 0
w0 = [x1, x2,z1,z2]
####----------------Function model passed to the ode solver--------------------
def f(t):
result = 0
for i in range(-25,25,1):
result = result + 1.0*DiracDelta(t-(i+1)*2*np.pi)
def vectorfield(t, w, omg1, omg2,Omg,Omegar, beta,Mu,K,g,m):
x1, x2,z1,z2= w
#result = 0
#for i in range(-10,10,1):
# result = result + 1.0*DiracDelta(t-(i+1)*np.pi)
field = [((g*np.sqrt(2*(x1)/(m*Omegar))*np.sin(x2))*(-0.5* abs(z2)**2 * np.cos(np.pi/3) + 0.5* np.conj(z1)* z2* np.sin(np.pi/3) + 0.5*
abs(z1)**2*np.cos(np.pi/3) + 0.5* np.conj(z2)* z1* np.sin(np.pi/3))-K*np.sin(x2)*f(t)),
((-g/np.sqrt((2*(x1)*m*Omegar))*np.cos(x2))*(-0.5* abs(z2)**2 * np.cos(np.pi/3) + 0.5* np.conj(z1)* z2* np.sin(np.pi/3) + 0.5*
abs(z1)**2*np.cos(np.pi/3) + 0.5* np.conj(z2)* z1* np.sin(np.pi/3))+x1),
((z1*(0.5*Omg+0.5*g*np.sqrt(2*x1/(m*Omegar))*np.cos(x2)*np.cos(np.pi/3))+0.5*z2*g*np.sqrt(2*x1/ (m*Omegar))*np.cos(x2)*np.sin(np.pi/3)) * -1j * (1/hcut) * (1/(cmath.sqrt(abs(z1)**2 + abs(z2)**2 )))),
(( z2*(-0.5*Omg+0.5*g*np.sqrt(2*x1/(m*Omegar))*np.cos(x2)*np.cos(np.pi/3))+ 0.5*g*np.sqrt(2*x1/(m*Omegar))*np.cos(x2)*z1*np.sin(np.pi/3)) * -1j * (1/hcut) * (1/(cmath.sqrt(abs(z1)**2 + abs(z2)**2 ))))]
#field1 = np.array(field, dtype='complex_')
#print(abs(z1)**2 + abs(z2)**2 )
print(z2)
return field
duration = 50
# time points
t = np.linspace(0, duration, 100)
abserr = 1.0e-10
relerr = 1.0e-6
#solution = odeint(vectorfield, w0, t, args=(p,))
solution = solve_ivp(vectorfield, [0, duration], w0,t_eval=t ,args=(p), atol=abserr,
rtol=relerr)
lw = 1
#'''
plot1 = plt.figure(1)
plt.style.use('seaborn-darkgrid')
plt.xlabel('time(t)')
plt.grid(True)
####----------------Plotting the oscillator dynamics---------------------------
plt.plot(t, solution.y[0,:], 'b', label='I', linewidth=lw)
plt.plot(t, solution.y[1,:], 'r', label='$\Theta$', linewidth=lw)
plt.plot(t, solution.y[2,:], 'g', label='x1(t)', linewidth=lw)
plt.plot(t, solution.y[3,:], 'orange', label='x2(t)', linewidth=lw)
plt.legend()
expect_1 = np.absolute(solution.y[2,:])**2 - np.absolute(solution.y[3,:])**2
plot2 = plt.figure(2)
plt.xlabel('time(t)')
#plt.plot(t, result, 'm', label='z_expect2', linewidth=lw)
plt.plot(t, solution.y[0,:], 'b', label='I', linewidth=lw)
plt.legend()
plt.show()
目前没有回答
相关问题 更多 >
编程相关推荐