python - 我想在 python 代码中使用 Dirac delta 作为时间函数来求解四个耦合微分方程。我没有做对

Question

我想在 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()

标签： pythonscipydifferential-equationsodeint