Is it possible to solve difference equation in MATLAB?

Can I solve this below difference equation using Matlab?

y(k)-3y(k-1)+2y(k-2)=2u(k-1)-2u(k-2)

u(k)=k when k>=0 and u(k)=0 when k<0

y(k)=0 when k<0

[John D'Errico]

A cute homework problem. But you would find it difficult to convince me it is not HW. Anyway, I'll show a complete solution in time, at least if anyone cares.

Answer 1

[John D'Errico]

It looks like this is a homework question, so I am reluctant to post a complete answer. But at the same time, it looks like you already want to use filter. Therefore it won't be doing your homework (because what matters is how I arrived at that form) if I suggest this to consider:

Y = @(k) -4 - 2*k + 4*2.^k;
Y((0:10)')

ans = 11×1

0 2 8 22 52 114 240 494 1004 2026

The question is then, how did I arrive at that conclusion?

As a hint, consider the solution of the homogeneous problem, thus

y(k) - 3*y(k-1) + 2*y(k-2) == 0

So what functional form will that 3 term difference universally kill off? (This solution is not that different from how one would generate the Binet solution for the Fibonacci number sequence.)

Once you have the solution to that problem, now look at u. What does it contribute? In the end, you would come to the same conclusion as I did about the form for y(k).

Can you use filter? Well, yes, you should be able to do that. First, compute the sequence u. It is known, and trivial to compute. Then call filter. (u will be x in the call to filter.) I'll try not to say any more than that, since again, that would be doing homework. Of course, a loop itself is even simpler to write, because that requires no thinking at all. But soooooo boring to use loops.

For those of you who want to see how I derived the analytical form above, or who want to see the filter call that can replicate that solution for a finite sequence of numbers, I'll return in a few days to this question, expanding my answer in some detail then. But not immediately, just in case this was homework.

[Paul]

@John D'Errico

Can you double check your result? It doesn't match the result of @Benjamin Thompson in his Answer, which matches the result I derived (not shown).

[Tasin Nusrat]

It's not a homework as it's written in the first place that we cannot tell here to solve a complete project. But yes it's a book problem where students usually solve it 70% by hand calculation using equation and rest by matlab code. But as per my understanding I should be able to solve such kind of problems completely by MATLAB using filter or filtic function. The only problem I am getting is I don't know how to feed the conditions. As in,

clc;

clear all;

close all;

n = 0:10;

a = [1 -3 2]; % right hand side of difference equation

b = [0 2 -2]; % right hand side of difference equation

ui=n; %for n>=0 but how to write the other condition when n<0

y = filter(b,a,ui)

stem(n,y)

But thanks a lot for your detailed explianation. I will be looking forward to see how you have done it with filter command, whenever you feel to post it. Thank you!

[John D'Errico]

Drat! You were so close to the correct use of filter too. I'll guess that gets you close enough to let me show how to use filter for this.

help filter

 FILTER One-dimensional digital filter.
    Y = FILTER(B,A,X) filters the data in vector X with the
    filter described by vectors A and B to create the filtered
    data Y.  The filter is a "Direct Form II Transposed"
    implementation of the standard difference equation:
 
    a(1)*y(n) = b(1)*x(n) + b(2)*x(n-1) + ... + b(nb+1)*x(n-nb)
                          - a(2)*y(n-1) - ... - a(na+1)*y(n-na)
 
    If a(1) is not equal to 1, FILTER normalizes the filter
    coefficients by a(1). 
 
    FILTER always operates along the first non-singleton dimension,
    namely dimension 1 for column vectors and non-trivial matrices,
    and dimension 2 for row vectors.
 
    [Y,Zf] = FILTER(B,A,X,Zi) gives access to initial and final
    conditions, Zi and Zf, of the delays.  Zi is a vector of length
    MAX(LENGTH(A),LENGTH(B))-1, or an array with the leading dimension 
    of size MAX(LENGTH(A),LENGTH(B))-1 and with remaining dimensions 
    matching those of X.
 
    FILTER(B,A,X,[],DIM) or FILTER(B,A,X,Zi,DIM) operates along the
    dimension DIM.
 
    Tip:  If you have the Signal Processing Toolbox, you can design a
    filter, D, using DESIGNFILT.  Then you can use Y = FILTER(D,X) to
    filter your data.
 
    See also FILTER2, FILTFILT, FILTIC, DESIGNFILT.
 
    Note: FILTFILT, FILTIC and DESIGNFILT are in the Signal Processing
    Toolbox.

    Documentation for filter
       doc filter

    Other functions named filter

       arima/filter            gjr/filter         statespace/filter
       codistributed/filter    gpuArray/filter    tall/filter
       dssm/filter             LagOp/filter       timeseries/filter
       egarch/filter           msVAR/filter       varm/filter
       fints/filter            regARIMA/filter    vecm/filter
       garch/filter            ssm/filter

Here, I think the trick you were missing is that x in filter is just u. And we know u. Say we want to find the first 10 values. u is known of course.

u = (1:10)'; % transposed to make it easier to read the result.
a = [1 -3 2]; % you had this correct.
b = [0 2 -2]; % this as well.
filter(b,a,u)

ans = 10×1

0 2 8 22 52 114 240 494 1004 2026

Filter implicity assumes zero for the values that fall off the end of its universe.

The analytical solution just uses the method of undetermined coefficients, where I first derived a general solution for the homogeneous problem, then find a particular solution given u. So it looks a bit like the solution of an ODE, which it should resemble.

[Tasin Nusrat]

Thanks a lot!

[Paul]

The output of filter for the input u1 matches the ouptut of the closed form solution for Y in the original answer. However, that closed form solution does not satisfy the original difference equation and initial conditions in the problem statement. The difference equation is:

y(k) = 3y(k-1) - 2y(k-2) + 2u(k-1) - 2u(k-2) ; y(k) = 0 for k < 0; u(k) = 0 for k < -0; u(k) = k for k >= 0.

So, as @Benjamin Thompson showed

y(0) = 3*y(-1) - 2*y(-2) + 2*u(-1) - 2*u(-2) = 0

y(1) = 3*y(0) - 2*y(-1) + 2*u(0) - 2*u(-1) = 3*0 - 2*0 + 2*0 - 2*0 = 0

But the closed form solution (and that filter() ouput) is

Y(1) = -4 - 2*1 + 4*2.^1 = -6 + 8 = 2

[John D'Errico]

Sorry. But no. Filter returns the correct sequence. It starts at y(1). This is consistent with what others have said.

Admittedly, I fell asleep when I computed the corresponding coefficients for y(k), using the sequence starting one spot off, and therefore my coefficients in y(K) were wrong. It should have read:

y = @(k) -2*k - 2 + 2*2.^k;
y((1:10)')

ans = 10×1

0 2 8 22 52 114 240 494 1004 2026

[Paul]

Now that we have the correct expression for the closed form solution, we can see that @Tasin Nusrat did use filter() correctly in this comment. Comparing the two solutions:

y = @(k) -2*k - 2 + 2*2.^k;
n = 0:10;
a = [1 -3 2];  % right hand side of difference equation
b = [0 2 -2];  % right hand side of difference equation
    ui=n; %for n>=0 but how to write the other condition when n<0
yfilt = filter(b,a,ui)

yfilt = 1×11

0 0 2 8 22 52 114 240 494 1004 2026

isequal(y(0:10),yfilt)

ans = logical
   1

Why do you think that @Tasin Nusrat's solution was only "close to the correct use of filter"?

Answer 2

[Paul]

At least three ways this problem can be solved in Matlab that have been discussed in this thread. I thought it might be helpful to put them all together.

Restatement of the problem:

y(k)-3y(k-1)+2y(k-2)=2u(k-1)-2u(k-2)

u(k)=k when k>=0 and u(k)=0 when k<0

y(k)=0 when k<0

Option 1: Because the initial conditions on the output are zero and the input is causal, we can use filter(), exactly like @Tasin Nusrat did to solve for the first 11 outputs of y

k = 0:10;
a = [1 -3 2];  % left hand side of difference equation
b = [0 2 -2];  % right hand side of difference equation
u = k; % for k >=0
yfilter = filter(b,a,u).'; % yfilter(1) corresponds to k = 0, etc.

Option 2: Use a loop, credit to @Benjamin Thompson. This code could be cleaned up, but I think it has the benefit of being clear. Note that we have to account for Matlab's 1-based indexing.

yloop = nan(11,1);
u = @(k) (k.*(k>=0));
for kk = 0:10
    uk1 = u(kk-1);
    uk2 = u(kk-2);
    if kk == 0
        yk1 = 0;
        yk2 = 0;
    elseif kk == 1
        yk1 = yloop(0 + 1);
        yk2 = 0;
    else
        yk1 = yloop(kk-1 + 1);
        yk2 = yloop(kk-2 + 1);
    end
    yloop(kk + 1) = 3*yk1 - 2*yk2 + 2*uk1 - 2*uk2;
end

Option 3: Closed form solution approach, suggested by @John D'Errico, but using z-transform.

syms k integer
syms z

% Transfer function from U to Y, by inspection
H(z) = (2/z - 2/z^2)/(1 - 3/z + 2/z^2);

% z-transform of U(z)
syms U(z)
U(z) = ztrans(k,k,z);

% Because the initial conditions on y are zero, and because u is causal, the z-transform of the output is
syms Y(z)
Y(z) = simplify(H(z)*U(z));

% closed form expression for y(k)
y(k) = simplify(iztrans(Y(z),z,k))

y(k) = 

% make sure y(k) = 0 for k < 0
y(k) = piecewise(k<0, 0, y(k));

% evaluate
yclosed = double(y(0:10)).';

Verify all three solutions are the same:

k = (0:10).';
table(k, yfilter , yloop , yclosed)

ans = 11×4 table
    k     yfilter    yloop    yclosed
    __    _______    _____    _______

     0        0         0         0  
     1        0         0         0  
     2        2         2         2  
     3        8         8         8  
     4       22        22        22  
     5       52        52        52  
     6      114       114       114  
     7      240       240       240  
     8      494       494       494  
     9     1004      1004      1004  
    10     2026      2026      2026  

We know that y = 0 for k < 0, so we can prepend yloop or yfilter with zeros to plot values on the negative axis if desired.

stem(-5:10,[zeros(1,5) yfilter.'])

Answer 3

[Benjamin Thompson]

Yes it is. Reorder terms so y(k) is on one side of the equation. Write a for loop in MATLAB to solve this numerically after you decide how many steps in k you want to solve for.

[Tasin Nusrat]

How can I define the conditions for u(k) and y(k)? And also can i do it using filter function?

[Benjamin Thompson]

You must solve it implicitly for values of k you are interested in. The sequence for u(k) is already defined per your original question. So lets say for 0 <= k < 10

y(0) = 0

y(1) = 3*y(0) - 2*y(-1) + 2*u(0) - 2*u(-1) = 0

y(2) = 3*y(1) - 2*y(0) + 2*u(1) - 2*u(0) = 2

y(3) = and repeat...

You can write a for loop in MATLAB or in many other programming languages to solve this for additional values of k.

[John D'Errico]

Actually, you don't NEED to go to that extent. As I say in my answer, you can indeed use filter, or you can find an explicit functional form, though that took me a little thought to compute the actual coefficients of such a form. If you want to see how to use filter, I'll come back and show that.