Calculate the first n terms of the recursion

2 views (last 30 days)
Calculate the first n terms of the following recursion
function xn_plus_one = state(n,t)
if n<1 || ~isscalar(n) || fix(n)~=n
error('n must be a non-negative integer');
else
xn_plus_one = zeros(1,length(t));
for jj = 1:length(t)
if t(jj) < 0 || t(jj) > 3
error('t must contain non-negative scalars less than or equal to 3');
else
if n==1
xn_plus_one(jj)=1;
else
xn_plus_one(jj) = integral(@(s) (2.*control(n-1,s)),0,t(jj)-1);
end
end
end
end
end
function C=control(n,v)
C=zeros(1,length(v));
for ii=1:length(v)
if v(ii)>=-1 && v(ii)<0
C(ii)=0;
elseif v(ii)>=0 && v(ii)<=2
C(ii)=integral(@(s) state(n,s),0,v(ii));
elseif v(ii)>2 && v(ii)<=3
C(ii)=0;
else
error('entries in v must vary from -1 to 3')
end
end
end
Calling the state function at n=4 and t=3 takes too much time to run. How can I improve it?

Accepted Answer

Morgan
Morgan on 7 Feb 2024
Edited: Morgan on 7 Feb 2024
Using the built-in MATLAB profiler, the biggest computational load is in the evaluations of integral(). Which is no surprise since it is notoriously slow inside for loops, let alone inside a double for-loop inside a recursive function, when there are much more efficient ways of integration.
Below, I've rewritten your code to exclude for loops and calculate everything as efficiently as I could come up with. Also, I shifted n down by one compared to your code to be more consistent with the equations you provided.
% INITIALIZE MATLAB
clear variables
close all
clc
% DEFINE INPUT ARGUMENTS
n = 1;
t = -2 : 0.01 : 5;
% CALL RECURSIVE STATE FUNCTION
x = state(n,t);
% PLOT X
figure(1);
plot(t,x,'-b','LineWidth',2);
% STATE FUNCTION
function x = state(n,t)
% CHECK INPUT ARGUMENTS
if n < 0 || ~isscalar(n) || fix(n)~=n
error('"n" must be a non-negative integer.');
end
% INITIALIZE X
N = length(t);
x(1,:) = ones(1,N);
% BASE CASE
if n == 0
return
end
% CALCULATE VALID INDICES OF t FOR x
idx1 = t >= 0 & t <= 3;
% ALL OTHER CASES
x(2:n+1,:) = zeros(n,N);
for i = 2 : n+1
% CALCULATE u
u = control(i-1,t-1); % <-- Edit: upper bound t-1 not t
% CALCULATE STATE FUNCTION
x(i,idx1) = 2*cumtrapz(t(idx1),u(idx1));
end
% RETURN FINAL ROW
x = x(n+1,:);
function u = control(n,t)
% ASSUME u IS ZERO
u = zeros(size(t));
% FIND VALID INDICES OF t FOR u
idx2 = t >= 0 & t <= 2;
% CALCULATE CONTROL FUNCTION
u(idx2) = cumtrapz(t(idx2), x(n,idx2));
end
end
Edit: My apologies, I forgot the upper bound on the integral of . It is highlighted in the comments above.

More Answers (0)

Categories

Find more on Mathematics in Help Center and File Exchange

Community Treasure Hunt

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

Start Hunting!