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

问题描述

我想在python代码中使用Dirac delta作为时间的函数来解决四个耦合 微分方程。在python代码中，我使用solve_ivp来求解耦合方程。
我不明白。我无法定义 Dirac delta 函数。请帮助我，如果 任何人都行。在python代码中，我使用solve_ivp来求解耦合方程。
我不明白。我无法定义 Dirac delta 函数。请帮助我，如果 任何人都可以。

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,m):
    
    x1,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],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$',solution.y[2,'g',label='x1(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.legend()

plt.show()

解决方法

暂无找到可以解决该程序问题的有效方法，小编努力寻找整理中！

如果你已经找到好的解决方法，欢迎将解决方案带上本链接一起发送给小编。

小编邮箱:dio#foxmail.com (将#修改为@）