0

I want to use Dirac delta as function of time in a python code to solve four coupled differential equation. In the python code i am using solve_ivp to solve the coupled equation.
I am not getting it right. I am not able to define the Dirac delta function. Please help me if anyone can. In the python code i am using solve_ivp to solve the coupled equation.
I am not getting it right. I am not able to define the Dirac delta function. Please help me if anyone can.

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()
user68207
  • 23
  • 4
  • 2
    Please make a minimal example that demonstrates what kind of effect you want to achieve. Also, consult your sources, the Dirac delta is a distribution, a functional on a function space. At most it can be interpreted as the mass density of a point mass. It never is an ordinary function. You will need to simulate the effect of delta by interrupting the integration and manually inserting the jump in the appropriate components. – Lutz Lehmann Mar 25 '21 at 18:12
  • Dear Sir, Yes, I exactly needs this what you said but I don't know how to do it I have Dirac delta in the code F(t) will denote as kick for the equation when (t= 2*n *Pi) .( You will need to simulate the effect of delta by interrupting the integration and manually inserting the jump in the appropriate components.) – user68207 Mar 25 '21 at 18:43
  • I see a call to `DiracDelta` in `f`, but don't see any kind of definition or import. – hpaulj Mar 25 '21 at 19:12
  • See https://stackoverflow.com/questions/63230568/numerical-solution-to-a-differential-equation-containing-a-dirac-delta, I have put two variants to deal with the delta in a new answer. The ODE there is minimal, but the deltas might not be in a place where you would need them. – Lutz Lehmann Mar 25 '21 at 19:20

0 Answers0