Converting Differential to Difference Equation - Frequency Domain Analysis

Calculus Level 5

Consider the following differential equation:

x ˙ = x + u \dot{x} = x + u

Where x x is a time-varying dynamic output quantity and u u is an input signal. The transfer function associated with this input-output relationship is G c ( s ) G_c(s) . This differential equation is converted to a discrete-time difference equation using three methods. The sampling time step is h = 0.1 h = 0.1 seconds.

Method 1:

x ( k + 1 ) = ( 1 + h ) x ( k ) + h u ( k ) x(k+1) = (1+h) x(k) + h u(k)

The transfer function associated with this input-output relationship is G 1 ( z ) G_1(z) .

Method 2:

x ( k + 1 ) = e h x ( k ) + ( e h − 1 ) u ( k ) x(k+1) = e^h x(k) + (e^h-1) u(k)

The transfer function associated with this input-output relationship is G 2 ( z ) G_2(z) .

Method 3:

The transfer function associated with this input-output relationship is G 3 ( z ) G_3(z) . G 3 ( z ) G_3(z) is obtained by replacing the Laplace-domain variable s s in G c ( s ) G_c(s) , by the following:

s = 2 h ( z − 1 z + 1 ) s = \frac{2}{h}\left(\frac{z-1}{z+1}\right)

It is to be noted that s = j ω s = j\omega and z = e j ω h z = e^{j \omega h} . The frequency ω \omega is varied from 0 0 to 20 20 rad/sec. By making these substitutions, the transfer functions become frequency dependent complex numbers. The phase of these complex numbers can be found using the relationship:

P = arctan ⁡ ( I m ( G ) R e ( G ) ) P = \arctan\left(\frac{Im(G)}{Re(G)}\right)

Let the phase associated with G c ( s ) G_c(s) , G 1 ( z ) G_1(z) , G 2 ( z ) G_2(z) , G 3 ( z ) G_3(z) be P c ( ω ) P_c(\omega) , P 1 ( ω ) P_1(\omega) , P 2 ( ω ) P_2(\omega) , P 3 ( ω ) P_3(\omega) respectively.

Compute the area bounded by the curve P i ( ω ) P_i(\omega) , and the lines ω = 0 \omega = 0 , ω = 20 \omega = 20 and P i = 0 P_i = 0 . Let this area be A i A_i . Here, i ∈ { 1 , 2 , 3 } i \in \{1,2,3\} . Also compute the similar area for the phase function P c ( ω ) P_c(\omega) . Let this area be A c A_c .

Let:

R i = ∣ A c − A i ∣ R_i = \mid A_c - A_i \mid

Thus, three values R 1 R_1 , R 2 R_2 and R 3 R_3 are obtained. Enter your answer as the minimum among these three values .

Notes/Hints:

  • Use the Laplace transform to find the transfer function for the differential equation.

  • Use the z-transform to obtain the transfer functions for the difference equations

  • j = − 1 j = \sqrt{-1}

Bonus:

Based on this analysis, which method is most suitable for converting a differential to a difference equation? This question focusses on phase. What comments can be made on magnitude?

Problem based on a suggestion by Steven Chase

Also try:

Converting Differential to Difference Equation


The answer is 0.170214.

Karan Chatrath by Karan Chatrath , 27, Netherlands

This section requires Javascript.
You are seeing this because something didn't load right. We suggest you, (a) try refreshing the page, (b) enabling javascript if it is disabled on your browser and, finally, (c) loading the non-javascript version of this page . We're sorry about the hassle.

1 solution

Steven Chase
Jul 26, 2019

I will use a slightly different convention for the difference equation, but it is equivalent.

Methods 1 and 2

Methods 1 1 and 2 2 have the same form:

x k = a x k − 1 + b u k − 1 a 1 = 1 + h b 1 = h a 2 = e h b 2 = e h − 1 x_k = a \, x_{k-1} + b \, u_{k-1} \\ a_1 = 1 + h \\ b_1 = h \\ a_2 = e^h \\ b_2 = e^h - 1

Taking the z-transform of both sides:

X = a z − 1 X + b z − 1 U X ( 1 − a z − 1 ) = b z − 1 U H = X U = b z − 1 1 − a z − 1 = b z − a z = e j ω h X = a \, z^{-1} X + b \, z^{-1} U \\ X (1 - a \, z^{-1}) = b \, z^{-1} U \\ H = \frac{X}{U} = \frac{b \, z^{-1}}{1 - a \, z^{-1}} = \frac{b }{z - a } \\ z = e^{j \omega \, h}

The transfer function H H is a complex number, with a polar representation consisting of a magnitude and a phase angle.

Methods 3 and "c"

Start with the differential equation, and take the Laplace transform

x ˙ = x + u s X = X + U X ( s − 1 ) = U H = X U = 1 s − 1 \dot{x} = x + u \\ s X = X + U \\ X (s-1) = U \\ H = \frac{X}{U} = \frac{1}{s-1}

Per the problem statement, there are two candidate values for s s : the classic representation j ω j \omega and another representation in terms of the complex number z z .

s c = j ω s 3 = 2 h ( z − 1 z + 1 ) s_c = j \omega \\ s_3 = \frac{2}{h} \Big( \frac{z - 1}{z + 1} \Big)

Code and Plots

The attached code plots the magnitude and phase responses for the different methods against the angular frequency. It also calculates the required areas, as per the problem requirement. Method "c" gives the classic Laplace result, and is taken as a reference quantity.

What stands out to me is that relative to Methods 1 and 2, Method 3 has a slightly inferior magnitude response, and a substantially superior phase response. 

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
 
import math
import cmath

h = 0.1

a1 = 1.0+h
b1 = h

a2 = math.exp(h)
b2 = math.exp(h) - 1.0

dw = 10.0**(-4.0)

##############################################################

A1 = 0.0
A2 = 0.0
A3 = 0.0
Ac = 0.0

w = 0.0

#print "w M1 M2 M3 Mc P1 P2 P3 Pc"

while w <= 20.0:

    z = complex(math.cos(w*h),math.sin(w*h))

    s3 = (2.0/h)*((z-1.0)/(z+1.0))
    sc = complex(0.0,w)

    H1 = b1/(z-a1)
    H2 = b2/(z-a2)
    H3 = 1.0/(s3-1.0)
    Hc = 1.0/(sc-1.0)

    M1 = abs(H1)
    M2 = abs(H2)
    M3 = abs(H3)
    Mc = abs(Hc)

    P1 = cmath.phase(H1)
    P2 = cmath.phase(H2)
    P3 = cmath.phase(H3)
    Pc = cmath.phase(Hc)

    dA1 = P1 * dw
    dA2 = P2 * dw
    dA3 = P3 * dw
    dAc = Pc * dw

    A1 = A1 + dA1
    A2 = A2 + dA2
    A3 = A3 + dA3
    Ac = Ac + dAc

    #print w,M1,M2,M3,Mc,P1,P2,P3,Pc

    w = w + dw

##############################################################

R1 = math.fabs(Ac - A1)
R2 = math.fabs(Ac - A2)
R3 = math.fabs(Ac - A3)

Q = min(R1,R2,R3)

print ""
print ""
print A1
print A2
print A3
print Ac
print ""
print Q

# Results

#-45.1062142111
#-45.2397915274
#-35.2420984481
#-35.4123132822

#0.170214834115

2 Helpful 1 Interesting 1 Brilliant 0 Confused

Very nice solution, thank you!

Even I made similar observations about the magnitude response. I would obviously choose method 3 over the other two, even though the magnitude response is a bit inferior. I also saw that for smaller time steps h h , the frequency response of the other two methods improves. But this happens at the expense of computational time.

  • Method 1 is based on the Derivative approximation as you highlighted.

  • Method 2 is based on the note I posted. The method is known as 'Zero-order hold'.

  • Method 3 is based on the trapezoidal method of numerical integration.

Karan Chatrath - 1 year, 10 months ago

Log in to reply

Indeed, trapezoidal rule is generally superior to Euler in terms of accuracy, so that makes sense. I like to use the "total vector error" metric, since it combines magnitude and phase together.

T V E = ∣ Expected vector − Actual vector ∣ ∣ Expected vector ∣ TVE = \frac{| \text{Expected vector} - \text{Actual vector}| }{| \text{Expected vector}|}

Steven Chase - 1 year, 10 months ago

Log in to reply

This is an interesting metric for comparison. Just to clarify, in this context, by the term vector, you mean the complex numbers themselves, right?

Karan Chatrath - 1 year, 10 months ago

Log in to reply

@Karan Chatrath – Yes, I'm calling it a vector, since I would visualize the complex numbers as vectors in the 2D plane

Steven Chase - 1 year, 10 months ago

It would also be interesting to look at this from a numerical stability perspective

Steven Chase - 1 year, 10 months ago

Log in to reply

Yes, I recall some of your problems on stability regions for different numerical solvers. This is more food for thought. Thanks

Karan Chatrath - 1 year, 10 months ago

0 pending reports

×

Problem Loading...

Note Loading...

Set Loading...