2

I want to compute this integral $\frac{1}{L}\int_{-\infty}^{t}H(t^{'})\exp(-\frac{R}{L}(t-t^{'}))dt^{'}$ using numpy.convolution, where $H(t)$ is heavside function. I am supposed to get this equals to $\exp(-\frac{R}{L}t)H(t)$ below is what I did, I changed the limitation of the integral into -inf to +inf by change of variable multiplying a different H(t) then I used this as my function to convolve with H(t)(the one inside the integral), but the output plot is definitely not a exp function, neither I could find any mistakes in my code, please help, any hint or suggestions will be appreciated!

import numpy as np
import matplotlib.pyplot as plt
R = 1e3 
L = 3. 

delta = 1
Nf = 100
Nw = 200
k = np.arange(0,Nw,delta)
dt = 0.1e-3 
tk = k*dt
Ng = Nf + Nw -2
n = np.arange(0,Nf+Nw-1,delta)
tn = n*dt

#define H
def H(n):
    H = np.ones(n)
    H[0] = 0.5
    return H

#build ftns that get convoluted
f = H(Nf)
w = np.exp((-R/L)*tk)*H(Nw)

#return the value of I
It = np.convolve(w,f)/L

#return the value of Voutput, b(t)
b = H(Ng+1) - R*It
plt.plot(tn,b,'o')
plt.show()
Improve this question
4
  • ehh the latex is not showing up and i cant add pictures Commented Jan 25, 2015 at 17:38
  • I assume you are trying to determine the resistor voltage response of an RL circuit to a step input? If you excuse the pun, your code looks a bit... convoluted. Commented Jan 26, 2015 at 17:52
  • Oh I tried to only take first half and then it works(and I multiplied It by dt which I was supposed to do since np.convolve assumes delta t =1 and this solved the units issue). however, I dont know why we only take the first half and works(Python suggests this may have something to do with boundary effects, but where do we decide where the cutoff is exactly?) Commented Jan 28, 2015 at 23:11
  • I don't understand what you mean by "only take the first half and works". Have a look at my code, and try and understand the graph. Commented Jan 28, 2015 at 23:20

2 Answers 2

Reset to default
4

The issue with your code is not so much programming as it is conceptual. Rewrite the convolution as Integral[HeavisideTheta[t-t']*Exp[-R/L * t'], -Inf, t] (that's Mathematica code) and upon inspection you find that H(t-t') is always 1 within the limits (except for at t'=t which is the integration limit... but that's not important). So in reality you're not actually performing a complete convolution... you're basically just taking half (or a third) of the convolution.

If you think of a convolution as inverting one sequence and then going one shift at the time and adding it all up (see http://en.wikipedia.org/wiki/Convolution#Derivations - Visual Explanation of Convolution) then what you want is the middle half... i.e. only when they're overlapping. You don't want the lead-in (4-th graph down: http://en.wikipedia.org/wiki/File:Convolution3.svg). You do want the lead-out.

Now the easiest way to fix your code is as such:

#build ftns that get convoluted
f = H(Nf)
w = np.exp((-R/L)*tk)*H(Nw)

#return the value of I
It = np.convolve(w,f)/L
max_ind = np.argmax(It)
print max_ind
It1 = It[max_ind:]

The lead-in is the only time when the convolution integral (technically sum in our case) increases... thus after the lead-in is finished the convolution integral follows Exp[-x]... so you tell python to only take values after the maximum is achieved.

#return the value of Voutput, b(t) works perfectly now!

Note: Since you need the lead-out you can't use np.convolve(a,b, mode = 'valid').

So It1 looks like: It1

b(t) using It1 looks like: enter image description here

There is no way you can ever get exp(-x) as the general form because the equation for b(t) is given by 1 - R*exp(-x)... It can't mathematically follow an exp(-x) form. At this point there are 3 things:

  1. The units don't really make sense... check them. The Heaviside function is 1 and R*It1 is about 10,000. I'm not sure this is an issue but just in case, the normalized curve looks as such: normalized

  2. You can get an exp(-x) form if you use b(t) = R*It1 - H(t)... the code for that is here (You might have to normalize depending on your needs):

    b = R*It1 - H(len(It1))
# print len(tn)
plt.plot(tn[:len(b)], b,'o')
plt.show()

And the plot looks like: exponential

  1. Your question might still not be resolved in which case you need to explain what exactly you think was wrong. With the info you've given me... b(t) can never have an Exp[-x] form unless the equation for b(t) is messed with. As it stands in your original code It1 follows Exp[-x] in form but b(t) cannot.
Improve this answer
Sign up to request clarification or add additional context in comments.

4 Comments

Sorry I don't quite get it. What do you mean by lead-in part? Isn't the code only counting the overlapping part? I used your method, but the plot b(t) is still not exp(-Rt/L) which is what it supposed to look like.(my plot seems like flipped now)
A convolution does not inherently only deal with cases where there's pure overlap. It also deals with cases such as (4-th graph down: en.wikipedia.org/wiki/File:Convolution3.svg). Note that the red curve is not completely within the blue curve... So the first few terms in the convolution deal with cases like that... when... two curves are not completely overlapping. I call this the lead-in but in essence its only the first few terms of the convolution.
The problem is about RL circuit response. a resistor R and an inductor L in series, we apply voltage H(t) across the pair in series(the input voltage), and need to find b(t) which is the voltage across the inductor alone. And the integral i mentioned in the beginning came from the derivation process, I already have the result b(t) = exp(-Rt/L)H(t) just by solving the differential equ, so if the plot from python is not like exp(-x) then i must have done the convolution in python wrong. Thanks I got what the lead-in means now, but shouldnt it count into convolution anyway?
I don't think the answer should be a decreasing exponential, just based on physical reasoning. If you apply a step voltage from 0 to 1V, then the steady state voltage across the resistor should just be that voltage, 1V, once transient effects die out. I think you might have an error in your integration, and it's more likely to be a 1-exp(-x) solution.
0

I think there's a bit of confusion here about convolution. We use convolution in the time domain to calculate the response of a linear system to an arbitrary input. To do this, we need to know the impulse response of the system. Be careful switching between continuous and discrete systems - see e.g. http://en.wikipedia.org/wiki/Impulse_invariance.

The (continuous) impulse response of your system (which I assume to be for the resistor voltage of an L-R circuit) I have defined for convenience as a function of time t: IR = lambda t: (R/L)*np.exp(-(R/L)*t) * H.

I have also assumed that your input is the Heaviside step function, which I've defined on the time interval [0, 1], for a timestep of 0.001 s.

When we convolve (discretely), we effectively flip one function around and slide it along the other one, multiplying corresponding values and then taking the sum. To use the continuous impulse response with a step function which actually comprises of a sequence of Dirac delta functions, we need to multiply the continuous impulse response by the time step dt, as described in the Wikipedia link above on impulse invariance. NB - setting H[0] = 0.5 is also important.

We can visualise this operation below. Any given red marker represents the response at a given time t, and is the "sum-product" of the green input and a flipped impulse response shifted to the right by t. I've tried to show this with a few grey impulse responses. enter image description here

The code to do the calculation is here.

import numpy as np
import matplotlib.pyplot as plt

R = 1e3     # Resistance
L = 3.      #Inductance
dt = 0.001  # Millisecond timestep

# Define interval 1 second long, interval dt
t = np.arange(0, 1, dt)

# Define step function
H = np.ones_like(t)
H[0] = 0.5              # Correction for impulse invariance (cf http://en.wikipedia.org/wiki/Impulse_invariance)

# RL circuit - resistor voltage impulse response (cf http://en.wikipedia.org/wiki/RL_circuit)
IR = lambda t: (R/L)*np.exp(-(R/L)*t) * H # Don't really need to multiply by H as w is zero for t < 0
# Response of resistor voltage
response = np.convolve(H, IR(t)*dt, 'full')

The extra code to make the plot is here:

# Define new, longer, time array for plotting response - must be same length as response, with step dt
tp = np.arange(len(response))* dt

plt.plot(0-t, IR(t), '-', label='Impulse response (flipped)')

for q in np.arange(0.01, 0.1, 0.01):
    plt.plot(q-t, IR(t), 'o-', markersize=3, color=str(10*q))
t  = np.arange(-1, 1, dt)
H = np.ones_like(t)
H[t<0] = 0.
plt.plot(t, H, 's', label='Unit step function')
plt.plot(tp, response, '-o', label='Response')
plt.tight_layout()
plt.grid()
plt.xlabel('Time (s)')
plt.ylabel('Voltage (V)')
plt.legend()
plt.show()

Finally, if you still have some confusion about convolution, I strongly recommend "Digital Signal Processing: A Practical Guide for Engineers and Scientists" by Steven W. Smith.

Improve this answer

2 Comments

thanks for answering first. you are right, its a RL circuit response problem. The input voltage H(t) is the voltage acrossthe pair R and L in series, and the step response b(t) (output voltage) is the voltage across the inductor alone.
Actually, having re-read your comment, you want the voltage across the inductor. If you check the Wikipedia page, the impulse response of the inductor voltage is a Dirac delta function minus the impulse response you used. Fix that and you will get the decaying exponential you expect.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.