# Digital Signal Processing

EXERCISES

Instructions:

1. (a) Modify the function ex with 2eqs to solve the IVP (L4.4) for 0 ≤ t ≤ 50 using the MATLAB
routine ode45. Call the new function LAB04ex1.
Let [t,Y] (note the upper case Y) be the output of ode45 and y and v the unknown functions.
Use the following commands to define the ODE:
function dYdt= f(t,Y)
y=Y(1); v=Y(2);
dYdt =

[v; sin(t)-7*v-5*y];
Plot y(t) and v(t) in the same window (do not use subplot), and the phase plot showing v
vsy in a separate window.
Add a legend to the first plot. (Note: to display v(t) = y’(t), use ’v(t)=y’’(t)’).
Add a grid. Use the command ylim([-1.2,1.2]) to adjust the y-limits for both plots.
Adjust the x-limits in the phase plot so as to reproduce the pictures in Figure L4g.
Include the M-file in your report. (b) For what (approximate) value(s) of t does y reach a local maximum in the window 0 t 50?
Check by reading the matrix Y and the vector t. Note that, because the M-file LAB04ex1.m
is a function file, all the variables are local and thus not available on the Command Window.
To read the matrix Y and the vector t, you need to modify the M-file by adding the line
[t Y].
(c) What seems to be the long term behavior of y?
(d) Modify the initial conditions to y(0) = 2, v(0) = 3 and run the file LAB04ex1.m with the
modified initial conditions. Based on the new graphs, determine whether the long term
behavior of the solution changes. Explain. Include the pictures with the modified initial

Nonlinear Problems
Nonlinear problems do not present any additional difficulty from an implementation point of view.

EXERCISES

2. (a) Consider the modified problem The ODE (L4.7) is very similar to (L4.4) except for the y2 term in the left-hand side. Because
of the factor y2 the ODE (L4.7) is nonlinear, while (L4.4) is linear. There is however very
little to change in the implementation of (L4.4) to solve (L4.7). In fact, the only thing that
needs to be modified is the ODE definition.
Modify the function defining the ODE in LAB04ex1.m. Call the revised file LAB04ex2.m. The
new function M-file should reproduce the pictures in Fig L4h.

(b) Compare the output of Figs L4g and L4h. Describe the changes in the behavior of the solutionin the short term.
(c) Compare the long term behavior of both problems (L4.4) and (L4.7), in particular the amplitude of oscillations.
(d) Modify LAB04ex2.m so that it solves (L4.7) using Euler’s method with N = 1000 in the
interval 0 ≤ t ≤ 50 (use the file euler.m from LAB 3 to implement Euler’s method; do
not delete the lines that implement ode45). Let [te,Ye] be the output of euler, and note that Ye is a matrix with two columns from which the Euler’s approximation to y(t) must be
extracted. Plot the approximation to the solution y(t) computed by ode45 (in black) and the
approximation computed by euler (in red) in the same window (you do not need to plot v(t)
nor the phase plot). Are the solutions identical? Comment. What happens if we increase the
value of N?

1. Solve numerically the IVP in the interval 0 ≤ t ≤ 50. Include the M-file in your report.
Is the behavior of the solution significantly different from that of the solution of (L4.7)?
Is MATLAB giving any warning message? Comment.

A Third-Order Problem

Consider the third-order IVP 1. (a) Write a function M-file that implements (L4.8) in the interval 0 ≤ t ≤ 50. Note that the
initial condition must now be in the form [y0,v0,w0] and the matrix Y, output of ode45,
has now three columns (from which y, v and wmust be extracted). On the same figure, plot
the three time series and, on a separate window, plot the phase plot using
figure(2); plot3(y,v,w,’k.-’);
grid on; view([-40,60])
xlabel(’y’); ylabel(’v=y’’’); zlabel(’w=y’’’’’); Do not forget to modify the function defining the ODE.
The output is shown in Fig. L4i. The limits in the vertical axis of the plot on the left were
deliberately set to the same ones as in Fig. L4h for comparison purposes using the MATLAB
command ylim([-1.2,1.2]).
Include the M-file in you report.
(b) Compare the output of Figs L4h and L4i. What do you notice?
(c) Differentiate the ODE in (L4.7) and compare to the ODE in (L4.8).
(d) Explain why the solution of (L4.7) also satisfies the initial conditions in (L4.8). Hint: Substitute t = 0 into (L4.7) and use the initial conditions for y and v.

Solution

euler.m

function [t,y] = euler(f,tspan,y0,N)

% Solves the IVP y’ = f(t,y), y(t0) = y0 in the time interval tspan = [t0,tf]

% using Euler’s method with N time steps.

% Input:

%  f     = name of inline function or function M-file that evaluates the ODE

%          (if not an inline function,  use: euler(@f,tspan,y0,N))

%          For a system, the f must be given as column vector.

%  tspan = [t0, tf] where t0 = initial time value and tf = final time value

%  y0    = initial value of the dependent variable. If solving a system,

%          initial conditions must be given as a vector.

%  N     = number of steps used.

% Output:

% t = vector of time values where the solution was computed

% y = vector of computed solution values.

m = length(y0);

t0 = tspan(1);

tf = tspan(2);

h = (tf-t0)/N;             % evaluate the time step size

t = linspace(t0,tf,N+1);    % create the vector of t values

y = zeros(m,N+1);   % allocate memory for the output y

y(:,1) = y0′;              % set initial condition

for n=1:N

y(:,n+1) = y(:,n) + h*f(t(n),y(:,n));   % implement Euler’s method

end

t = t’; y = y’;    % change t and y from row to column vectors

end

LAB04ex1.m

function LAB04ex1

t0 = 0; tf = 50; y0 = [-0.5;0.5];  % parameters for ode45()

[t,y] = ode45(@f,[t0,tf],y0);  % call of ode45

u1 = y(:,1); u2 = y(:,2);  % y in output has 2 columns corresponding to u1 and u2

figure(1);

plot(t,u1,’b-‘,t,u2,’r-‘); % plot the 2 outputs y(t) and y'(t)

grid on; % add a grid

legend(‘y(t)’,’v(t) = y”(t)’); % place a legend on the plot

ylim([-1.2 1.2]); % limit the y-axis

figure(2)

plot(u1,u2);

xlabel(‘u_1’); ylabel(‘u_2’); % plot the phase plot

xlabel(‘y’);ylabel(‘v = y”’); % plot the phase plot

grid on;

xlim([-0.5 1]);

ylim([-1.2 1.2]);

end

%———————————————————————-

functiondYdt= f(t,Y)

y=Y(1); v=Y(2);

dYdt =[v; sin(t)-7*v-5*y];

end

LAB04ex2.m

function LAB04ex2

t0 = 0; tf = 50; y0 = [-0.5;0.5]; % parameters for ode45()

[t,y] = ode45(@f,[t0,tf],y0); % call of ode45

u1 = y(:,1); u2 = y(:,2); % y in output has 2 columns corresponding to u1 and u2

figure(1); % open a new figure windoe

plot(t,u1,’b-‘,t,u2,’r-‘); % plot the 2 outputs y(t) and y'(t)

grid on; % add a grid

legend(‘y(t)’,’v(t) = y”(t)’); % place a legend on the plot

ylim([-1.2 1.2]); % limit the y-axis

figure(2); % open a new window

plot(u1,u2,’b-‘);axis square;xlabel(‘y’);ylabel(‘v = y”’); % plot the phase plot

grid on;

xlim([-0.5 1]);

ylim([-1.2 1.2]);

end

%———————————————————————-

functiondYdt = f(t,Y)

y = Y(1); v = Y(2);

dYdt = [v ; sin(t)-7*v*(y.^2)-5*y];

end

LAB04ex3.m

function LAB04ex3

t0 = 0; tf = 50; y0 = [-0.5;0.5]; % parameters for ode45()

[t,y] = ode45(@f,[t0,tf],y0); % call of ode45

u1 = y(:,1); u2 = y(:,2); % y in output has 2 columns corresponding to u1 and u2

figure(1); % open a new figure windoe

plot(t,u1,’b-‘,t,u2,’r-‘); % plot the 2 outputs y(t) and y'(t)

grid on; % add a grid

legend(‘y(t)’,’v(t) = y”(t)’); % place a legend on the plot

ylim([-1.2 1.2]); % limit the y-axis

figure(2); % open a new window

plot(u1,u2,’k-‘);axis square;xlabel(‘y’);ylabel(‘v = y”’); % plot the phase plot

grid on;

xlim([-0.5 1]);

ylim([-1.2 1.2]);

end

%———————————————————————-

functiondYdt = f(t,Y)

y = Y(1); v = Y(2);

dYdt = [v ; sin(t)-7*v*y-5*y];

end

LAB04ex4.m

function LAB04ex4

t0 = 0; tf = 50; y0 = [-0.5;0.5;1.625]; % parameters for ode45()

[t,Y] = ode45(@f,[t0,tf],y0); % call of ode45

y = Y(:,1); v = Y(:,2); w = Y(:,3);

figure(1); % open a new figure windoe

plot(t,y,’b-‘,t,v,’r-‘,t,w,’m-‘); % plot the 3 outputs y,y’, and y”

grid on; % add a grid

legend(‘y’,’v = y”’,’w = y””’); % place a legend on the plot

ylim([-1.2 1.2]); % limit the y-axis

figure(2); % open a new window

plot3(y,v,w,’k.-‘);

grid on; view([-40,60]);

xlabel(‘y’); ylabel(‘v = y”’); zlabel(‘w = y””’);

zlim([-4,2]); % limit the z-axis

xlim([-0.6 0.6]); % limit the y-axis

end

%———————————————————————-

functiondYdt = f(t,Y)

y = Y(1); v = Y(2); w = Y(3);

dYdt = [v ; w;cos(t)-7*(y.^2)*w-14*(v.^2)*y-5*v];

end

format compact

% MAT 275 MATLAB Assigment # 4

% Exercise 1

function LAB04ex1

t0 = 0; tf = 50; y0 = [-0.5;0.5];  % parameters for ode45()

[t,y] = ode45(@f,[t0,tf],y0);  % call of ode45

u1 = y(:,1); u2 = y(:,2);  % y in output has 2 columns corresponding to u1 and u2

figure(1);

plot(t,u1,’b-‘,t,u2,’r-‘); % plot the 2 outputs y(t) and y'(t)

legend(‘y(t)’,’v(t) = y”(t)’); % place a legend on the plot

ylim([-1.2 1.2]); % limit the y-axis

figure(2)

plot(u1,u2);

xlabel(‘u_1’); ylabel(‘u_2’); % plot the phase plot

xlabel(‘y’);ylabel(‘v = y”’); % plot the phase plot

gridon;

xlim([-0.5 1]);

ylim([-1.2 1.2]);

end

%———————————————————————-

functiondYdt= f(t,Y)

y=Y(1); v=Y(2);

dYdt =[v; sin(t)-7*v-5*y];

end

LAB04ex1 % call the LAB04ex1 function to solve an ODE  (b) The output y(t) reaches a local maximum at 15.1567 seconds and it’s 0.124.

(c) In the long term, the output y(t) behaves as a sinusoidal signal.

(d)

function LAB04ex1

t0 = 0; tf = 40; y0 = [2;3]; % parameters for ode45()

[t,Y] = ode45(@f,[t0,tf],y0); % call of ode45

u1 = Y(:,1); u2 = Y(:,2); % y in output has 2 columns corresponding to u1 and u2

figure(1); % open a new figure windoe

plot(t,u1,’b-‘,t,u2,’r-‘); % plot the 2 outputs y(t) and y'(t)

legend(‘y(t)’,’v(t) = y”(t)’); % place a legend on the plot

ylim([-1.5 1.5]); % limit the y-axis

figure(2); % open a new window

plot(u1,u2,’k-‘);axis square;xlabel(‘y’);ylabel(‘v = y”’); % plot the phase plot

gridon;

xlim([-1 1]);

ylim([-1.5 1.5]);

end

%———————————————————————-

functiondYdt = f(t,Y)

y = Y(1); v = Y(2);

dYdt = [v ; cos(t)-4*v-3*y];

end

LAB04ex1 % call the LAB04ex1 function to solve an ODE  It can be observed from the figures above that the long-term behavior of the solution doesn’t change when the initial conditions change. The solution will still behave as a sine wave in the long-term.

% Exercise 2

function LAB04ex2

t0 = 0; tf = 50; y0 = [-0.5;0.5]; % parameters for ode45()

[t,y] = ode45(@f,[t0,tf],y0); % call of ode45

u1 = y(:,1); u2 = y(:,2); % y in output has 2 columns corresponding to u1 and u2

figure(1); % open a new figure windoe

plot(t,u1,’b-‘,t,u2,’r-‘); % plot the 2 outputs y(t) and y'(t)

legend(‘y(t)’,’v(t) = y”(t)’); % place a legend on the plot

ylim([-1.2 1.2]); % limit the y-axis

figure(2); % open a new window

plot(u1,u2,’b-‘);axis square;xlabel(‘y’);ylabel(‘v = y”’); % plot the phase plot

gridon;

xlim([-0.5 1]);

ylim([-1.2 1.2]);

end

%———————————————————————-

functiondYdt = f(t,Y)

y = Y(1); v = Y(2);

dYdt =

[v ; sin(t)-7*v*(y.^2)-5*y]

;

end

LAB04ex2 % call the LAB04ex2 function to solve an ODE  (b) From the graphs obtained, it can be concluded that the solution y presents overshoot twice before finally reaching a steady-state after 26 seconds approximately. On the other hand, for the ODE solved in Exercise 1a, with the same initial conditions, the solution y presented a negative overshoot and finally reached to a steady-state after 15 seconds approximately. The difference that was observed for this particular nonlinear ODE is due to the term y2 that appears.

(c) Long-term, the solution of this ODE and the previous one behave as a sine wave but for the current ODE the amplitude is at 0.25 approximately whereas for the previous ODE the amplitude is at 0.125 approximately. This can be interpreted that for the nonlinear ODE the amplitude of the solution in the long-term has doubled.

(d)

function LAB04ex2

t0 = 0; tf = 50; y0 = [-0.5;0.5]; % parameters for ode45()

[t,y] = ode45(@f,[t0,tf],y0); % call of ode45

[te,Ye] = euler(@f,[t0,tf],y0,1000); % call of euler() function

u1 = y(:,1); % solution by ode45

u3 = Ye(:,1); % solution by euler

figure(1); % open a new figure windoe

plot(t,u1,’k-‘,te,u3,’r-‘); % plot y(t) from ode45() and euler()

legend(‘ode45()’,’euler()’); % place a legend on the plot

ylim([-1.2 1.2]); % limit the y-axis

end

%———————————————————————-

functiondYdt = f(t,Y)

y = Y(1); v = Y(2);

dYdt =[v ; sin(t)-7*v*(y.^2)-5*y];

end

LAB04ex2 % call the LAB04ex2 function to solve an ODE using ode45() and euler() It can be observed from the figure above that the solution of the ODE using the ode45() function and the euler() function are almost but not entirely identical. Increasing the value of N to 5000 will result in the following figure. From this graph it can be observed that both solutions produced are now the same. Therefore, increasing the value of N, the solution returned by the Euler method will be the same as the one produced by ode45().

% Exercise 3

function LAB04ex3

t0 = 0; tf = 50; y0 = [-0.5;0.5]; % parameters for ode45()

[t,y] = ode45(@f,[t0,tf],y0); % call of ode45

u1 = y(:,1); u2 = y(:,2); % y in output has 2 columns corresponding to u1 and u2

figure(1); % open a new figure windoe

plot(t,u1,’b-‘,t,u2,’r-‘); % plot the 2 outputs y(t) and y'(t)

legend(‘y(t)’,’v(t) = y”(t)’); % place a legend on the plot

ylim([-1.2 1.2]); % limit the y-axis

figure(2); % open a new window

plot(u1,u2,’k-‘);axis square;xlabel(‘y’);ylabel(‘v = y”’); % plot the phase plot

gridon;

xlim([-0.5 1]);

ylim([-1.2 1.2]);

end

%———————————————————————-

functiondYdt = f(t,Y)

y = Y(1); v = Y(2);

dYdt =

[v ; sin(t)-7*v*y-5*y]

;

end

LAB04ex3 % call the LAB04ex3 function to solve an ODE  From the plot of the solution (y(t) and y΄(t)) it can be observed that the solution is evaluated until t = 5.840631 and it differs significantly than the solution of the previous ODE. This means that the ODE is stiff and that there might be a vertical asymptote near that time. The message returned in MatlabCommandWindow is the following:

Warning: Failure at t=5.840631e+00.  Unable to meet integration tolerances without reducing the step size below thesmallest value allowed (1.421085e-14) at time t.

> In ode45 (line 308)

In LAB04ex3 (line 3)

% Exercise 4

function LAB04ex4

t0 = 0; tf = 50; y0 = [-0.5;0.5;1.625]; % parameters for ode45()

[t,Y] = ode45(@f,[t0,tf],y0); % call of ode45

y = Y(:,1); v = Y(:,2); w = Y(:,3);

figure(1); % open a new figure windoe

plot(t,y,’b-‘,t,v,’r-‘,t,w,’m-‘); % plot the 3 outputs y,y’, and y”

legend(‘y’,’v = y”’,’w = y””’); % place a legend on the plot

ylim([-1.2 1.2]); % limit the y-axis

figure(2); % open a new window

plot3(y,v,w,’k.-‘);

gridon; view([-40,60]);

xlabel(‘y’); ylabel(‘v = y”’); zlabel(‘w = y””’);

zlim([-4,2]); % limit the z-axis

xlim([-0.6 0.6]); % limit the y-axis

end

%———————————————————————-

functiondYdt = f(t,Y)

y = Y(1); v = Y(2); w = Y(3);

dYdt =

[v ; w;cos(t)-7*(y.^2)*w-14*(v.^2)*y-5*v]

;

end

LAB04ex4 % call the LAB04ex4 function to solve an ODE  (b) Comparing the solutions from this part with the one obtained in Exercise 2 graphically, that both of the solutions are the same.

(c) Let’s take the derivative of the ODE as defined in exercise 2 as below.    (d) Let’s solve the ODE of exercise 2 for as follows: For t = 0, the value of the previous relation is equal to:   This value is the initial condition of the ODE . Therefore, the solution of the ODE defined in exercise 2 satisfies the initial conditions of the ODE defined in this part.