Azzera filtri
Azzera filtri

Is it possible to solve difference equation in MATLAB?

68 visualizzazioni (ultimi 30 giorni)
Can I solve this below difference equation using Matlab?
u(k)=k when k>=0 and u(k)=0 when k<0
y(k)=0 when k<0
  1 Commento
John D'Errico
John D'Errico il 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.

Accedi per commentare.

Risposta accettata

John D'Errico
John D'Errico il 8 Feb 2022
Modificato: John D'Errico il 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;
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 Commenti
John D'Errico
John D'Errico il 8 Feb 2022
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;
ans = 10×1
0 2 8 22 52 114 240 494 1004 2026
Paul il 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
ans = logical
Why do you think that @Tasin Nusrat's solution was only "close to the correct use of filter"?

Accedi per commentare.

Più risposte (2)

Paul il 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
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 il 7 Feb 2022
Modificato: Benjamin Thompson il 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 Commenti
Benjamin Thompson
Benjamin Thompson il 8 Feb 2022
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
John D'Errico il 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.

Accedi per commentare.

Community Treasure Hunt

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

Start Hunting!

Translated by