Is it possible to solve difference equation in MATLAB?

133 views (last 30 days)
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
John D'Errico on 8 Feb 2022
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.

Sign in to comment.

Accepted Answer

John D'Errico
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)')
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.
  7 Comments
Paul
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)
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"?

Sign in to comment.

More Answers (2)

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

Benjamin Thompson
Benjamin Thompson on 7 Feb 2022
Edited: Benjamin Thompson on 7 Feb 2022
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.
  3 Comments
John D'Errico
John D'Errico on 8 Feb 2022
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.

Sign in to comment.

Community Treasure Hunt

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

Start Hunting!

Translated by