Is it possible to solve difference equation in MATLAB?

78 次查看（过去 30 天）
Tasin Nusrat
Tasin Nusrat 2022-2-7
评论： Paul 2022-2-9
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 个评论
John D'Errico
John D'Errico 2022-2-8
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.

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

John D'Errico
John D'Errico 2022-2-8
编辑：John D'Errico 2022-2-8
Ran in:
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 个评论
显示 5更早的评论
John D'Errico
John D'Errico 2022-2-8
Ran in:
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
Paul 2022-2-9
Ran in:
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"?

请先登录，再进行评论。

更多回答（2 个）

Paul
Paul 2022-2-8
Ran in:
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 2022-2-7
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 个评论
显示 1更早的评论
Benjamin Thompson
Benjamin Thompson 2022-2-8
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 2022-2-8
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.

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Calculus 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by