Calculate the first n terms of the recursion

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?

Answer 1

[Morgan]

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.