# Is it possible to solve difference equation in MATLAB?

133 views (last 30 days)

Show older comments

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

##### 1 Comment

John D'Errico
on 8 Feb 2022

### Accepted Answer

John D'Errico
on 8 Feb 2022

Edited: John D'Errico
on 8 Feb 2022

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)')

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.

##### 7 Comments

Paul
on 9 Feb 2022

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)

isequal(y(0:10),yfilt)

### More Answers (2)

Paul
on 8 Feb 2022

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

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))

% 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)

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.'])

##### 0 Comments

Benjamin Thompson
on 7 Feb 2022

Edited: Benjamin Thompson
on 7 Feb 2022

##### 3 Comments

John D'Errico
on 8 Feb 2022

### See Also

### Categories

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!