"Juan Miguel" wrote in message <ki7eia$2h7$1@newscl01ah.mathworks.com>...
> suppose I have a positive real number T, and I a have a vector of parameters "a" with N elements, and a function f(k,t) on t defined recursively in k, for k<N as follows
>
> (the following is NOT matlab code, it is just pseudo code to illustrate what I need to do, since i am not sure what Matlab data structure to use for the function)
>
> f(0,t) = 1 a constant for all t
>
> f(1,t) = @(t) integral( @(s) exp(a(1)*s),t,T);
> f(2,t) = @(t) integral( @(s) f(1,s)*exp(a(2)*s),t,T);
> .
> .
> .
> f(k,t) =@(t) integral( @(s) f(k1,s)*exp(a(k1)*s),t,T);
> ........
       
It isn't necessary to do iteration (recursion) to find your f(n,t) for a given n. It can be expressed directly as a linear combination of n+1 exponential functions. The following matlab code evaluates numerically the necessary parameters for such an expression.
     
clear
% Let a be an nelement column vector
a = .05+.1*rand(10,1); % < Choose your 'a' vector
% Use a to compute the scalar b and n+1element row vectors c and d.
n = size(a,1);
p = tril(reshape(1:n*(n+1),n,n+1));
p = p(p>0);
q = (n+1)*p(n*(n+1)+1)*floor(p/(n+1));
M = cumsum(tril(repmat(a,1,n+1)));
c = M(n,:);
b = c(1);
M(q) = M(p);
d = prod(M).*((1).^(0:n));
% Let T be a scalar and t an melement vector
T = rand;
m = 100;
t = Tlinspace(1,2,m);
% Compute the corresponding values of f(n,t) and plot it.
fn = zeros(size(t));
for ix = 1:m
fn(ix) = exp(b*T)*sum(exp(c*(Tt(ix)))./d);
end
plot(t,fn,'y') % < Plot of f(n,t) versus t
     
For example the corresponding expression for f(3,t) would be this:
f(3,t) = exp((a1+a2+a3)*T) * ...
( +exp((a1+a2+a3)*(Tt))/a1/(a1+a2)/(a1+a2+a3) ...
exp((a2+a3)*(Tt))/a1/a2/(a2+a3) ...
+exp(a3*(Tt))/a2/(a1+a2)/a3 ...
1/a3/(a2+a3)/(a1+a2+a3) )
where a1 = a(1), a2 = a(2), and a3 = a(3).
It is only fair to warn you that with a large value for n, even with double precision floating point accuracy, these numerical computations may encounter serious relative round off errors for values of t near T since it involves a linear combination of sizable exponential quantities whose sum should cancel out to nearly zero. The only remedy I can see for that is to use higher precision packages from the file exchange to do the equivalent of the above algorithm with the appropriate accuracy. If you try to carry this out using the symbolic toolbox, the expressions would have an order O(n^3) characterlength which means that if n is in the hundreds, the expression would have a characterlength in the millions!
One other warning for this algorithm is that no sequence of consecutive values from vector 'a' should have a sum of zero since that would yield an indeterminate ratio at some point which would result in a NaN value. That is, in the above example with n = 3, none of the values a1, a2, a3, a1+a2, a2+a3, and a1+a2+a3 can be zero.
Roger Stafford
