Declaring a functional recursive sequence in Matlab

Question

I'd like to declare first of all, that I'm a mathematician. This might be a stupid stupid question; but I've gone through all the matlab tutorials--they've gotten me nowhere. I imagine I could code this in C (it'd be exhausting); but I need matlab for this particular function. And I don't get exactly how to do it.

Here is the pasted Matlab code of where I'm running into trouble:

function y = TAU(z,n)
   y=0;
   for i =[1,n]
       y(z) = log(beta(z+1,i) + y(z+1)) - beta(z,i);
   end
end

(beta is an arbitrary "float" to "float" function with an index i.)

I'm having trouble declaring y as a function, in which we call the function at a different argument. I want to define y_n(z) with something something y_{n-1}(z+1). This is all done in a recursive process to create the function. I really feel like I'm missing something stupid.

As a default function it assigns y to be an array (or whatever you call the default index assignment). But I don't want an array. I want y to be assigned as a "function" class (i.e. takes "float" to "float"). And then I'm defining a sequence of y_n : "float" to "float". So that z to z+1 is a map on "float" to "float".

I don't know if I'm asking too much of matlab...

Help a poor mathematician who hasn't coded since the glory days of X-box mods.

...Please don't tell me I have to go back to Pari-GP/C drawing boards over something so stupid.

Please help!

EDIT: At rahnema1 & mimocha's request, I'll describe the math, and of what I am trying to do with my program. I can't see how to implement latex in here. So I'll write the latex code in a generator and upload a picture. I'm not so sure if there even is a work around to what I want to do.

As to the expected output. We'd want,

beta(z+1,i) + TAU(z+1,i) = exp(beta(z,i) + TAU(z,i+1))

And we want to grow i to a fixed value n. Again, I haven't programmed in forever, so I apologize if I'm speaking a little nonsensically.

EDIT2:

So, as @rahnema1 suggests; I should produce a reproducible example. In order to do this, I'll write the code for my beta function. It's surprisingly simple. This is for the case where the "multiplier" variable is set to log(2); but you don't need to worry about any of that.

function f = beta(z,n)
    f=0;
    for i = 0:n-1
        f = exp(f)/(1+exp(log(2)*(n-i-z)));
    end
end

This will work fine for z a float no greater than 4. Once you make z larger it'll start to overflow. So for example, if you put in,

beta(2,100)

 1.4242

beta(3,100)

 3.3235

beta(3,100) - exp(beta(2,100))/(1/4+1)

0

The significance of the 100, is simply how many iterations we perform; it converges fast so even setting this to 15 or so will still produce the same numerical accuracy. Now, the expected output I want for TAU is pretty straight forward,

TAU(z,1) = log(beta(z+1,1)) - beta(z,1)
TAU(z,2) = log(beta(z+1,2) + TAU(z+1,1)) - beta(z,2)
TAU(z,3) = log(beta(z+1,3) + TAU(z+1,2)) - beta(z,3)
...
TAU(z,n) = log(beta(z+1,n) + TAU(z+1,n-1)) -beta(z,n)

I hope this helps. I feel like there should be an easy way to program this sequence, and I must be missing something obvious; but maybe it's just not possible in Matlab.

At mimocha's suggestion, I'll look into tail-end recursion. I hope to god I don't have to go back to Pari-gp; but it looks like I may have to. Not looking forward to doing a deep dive on that language, lol.

Thanks, again!

Answer 1

Is this what you are looking for?

function out = tau(z,n)
    % Ends recursion when n == 1
    if n == 1
       out = log(beta(z+1,1)) - beta(z,1);
       return
    end

    out = log(beta(z+1,n) + tau(z+1,n-1)) - beta(z,n);
end

function f = beta(z,n)
    f = 0;
    for i = 0:n-1
        f = exp(f) / (1 + exp(log(2)*(n-i-z)));
    end
end

This is basically your code from the most recent edit, but I've added a simple catch in the tau function. I tried running your code and noticed that n gets decremented infinitely (no exit condition).

With the modification, the code runs successfully on my laptop for smaller integer values of n, where 1e5 > n >= 1; and for floating values of z, real and complex. So the code will unfortunately break for floating values of n, since I don't know what values to return for, say, tau(1,0) or tau(1,0.9). This should easily be fixable if you know the math though.

However, many of the values I get are NaNs or Infs. So I'm not sure if your original problem was Out of memory error (infinite recursion), or values blowing up to infinity / NaN (numerical stability issue).

Here is a quick 100x100 grid calculation I made with this code.

Then I tested on negative values of z, and found the imaginary part of the output to looks kinda cool.

Not to mention I'm slightly geeking out over the fact that pi is showing up in the imaginary part as well :)

tau(-0.3,2) == -1.45179335740446147085 +3.14159265358979311600i