Is it possible to solve difference equation in MATLAB?
133 views (last 30 days)
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;
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.
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:
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;
yk1 = yloop(kk-1 + 1);
yk2 = yloop(kk-2 + 1);
yloop(kk + 1) = 3*yk1 - 2*yk2 + 2*uk1 - 2*uk2;
Option 3: Closed form solution approach, suggested by @John D'Errico, but using z-transform.
syms k integer
% 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)
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
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));
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.