How to define Recurrence Relations in IGOR

Hi,

I don't have much experience with programming, but I have recurrence relations for which there is no closed-form. I want to define them in IGOR Pro. I used simple if..else to define them but, when I try to get an output I get the following error:

"out of stack space (too much macro or function recursion)". It just works for P0(0), P0(1), Q0(0) and Q0(1). What would be the way to go?

Please help. 

Function P0(n)

variable n

variable a,q
if(n==0)
return 0

elseif(n==1)
return 1

elseif(n>=2)
return (2*P0(n-1)-2*a*P0(n)-q*Q0(n))/(n+1)

endif

End

Function Q0(n)
variable n

variable a,q
if(n==0)
return 0

elseif(n==1)
return 0

elseif(n>=2)
return (2*Q0(n-1)-2*a*Q0(n)-q*P0(n))/(n)

endif

End

 

In your example above I can't see where the variables a and q are assigned any value?

Note also, q has a special meaning in Igor - it is the column index of a destination wave when used in a multidimensional wave assignment. See DisplayHelpTopic("Waveform Arithmetic and Assignment").

Both of your functions are (in the case of n >= 2) functions of themselves. i.e. the function calls itself *with the same parameter value*. This leads to an infinite nest of function calls. For example, calling P0(2) will (in the last return statement) calls a new instance of the function P0(2), which in turn calls a new instance of P0(2), and so on.

I have no idea whether the following is valid, but...

if you take the function assignment P0(n) = (2*P0(n-1)-2*a*P0(n)-q*Q0(n))/(n+1),

you can rearrange this to something like:

P0(n) = (2*P0(n-1) - q*Q0(n)) / (n + 1 + 2*a)

Does this help?

Regards,

Kurt

OP approach is less than ideal in that there are limits to recurrence (as OP found) and it is very inefficient because one is unnecessarily calculating previous values multiple times.  One approach is to estimate the largest parameter N that one wants to support and then create a couple of waves say:

Make/n=(N) pWave=nan,qWave=nan

From here define your functions similar to:

Function P0(n,pWave)
    variable n
    Wave pWave

    variable a,q,res

    if(numType(pWave[n])!=2)
        return pWave[n]
    endif
   
    if(n==0)
        pWave[0]=0
        return 0
    elseif(n==1)
        pWave[1]=1
        return 1
    elseif(n>=2)
        res=(2*P0(n-1,pWave)-2*a*P0(n,pWave)-q*Q0(n,pWave))/(n+1)
        pWave[n]=res
        return res
    endif

End

I hope this helps,

 

A.G.