Implement recursive definition and closed form in octave

Question

I want to find $I_n := \int_{0}^{1} x^n e^x dx$ with numerical methods.

I have found the recurrence relation $I_n = e - n \cdot I_{n - 1}$ and the closed form \begin{align} \tag{$\star$} \label{eq:closedform} I_n = \sum_{k = 0}^{n} (-1)^{k} \frac{n! e}{(n - k)!} - (-1)^n n!. \end{align} Now I want to implement two octave methods

1. I=intrek(n,I_0) with calculates $I_n$ with the recurrence relation
2. I=intexp(n), which calculates $I_n$ with the closed form \eqref{eq:closedform}.

I know python but am finding octave syntax very unintuitive, but I know this should be to difficult since it won't involve more complex structures that i.e while-loops.

Here's what I tried so far.

function I=intrek(n,I0)
  while (i < n + 1)
    x = exp(1) - i*I0;
    i++;
  endwhile
  x
endfunction

function I = intexp(n)
  x = sum(arrayfun(@(k) (-1)^k * (n! / (n - k!)) * exp(1), 1:5));
  y = x - (-1)^n * n!
endfunction

Disclaimer This are homework from previous courses I am trying to solve for exam preparation purposes.

I then have to "test" my function with the following code:

for k=1:5
    n=10*k;
    I=0;
    h=1/(2*10000);
    for i=0:10000-1
        x=i*2*h;
        I=I+x^n*exp(x);
        x=x+h;
        I=I+4*x^n*exp(x);
        x=x+h;
        I=I+x^n*exp(x);
    end
    corr(k)=h/3*I;
    rek(k)=intrek(n,exp(1)-1);
    expl(k)=intexp(n);
end
disp(' ')
disp(['n        ','exact                 ', 'recursive               ', 'explicit'])
disp(' ')
res=[(10:10:50)',corr',rek',expl'];
disp(num2str(res))
!rm results.txt
!rm output.txt
!rm results.ps
save -ASCII 'results.txt' res
!more results.txt intrek.m intexp.m>output.txt
!mpage -2A  output.txt> results.ps
close all
figure
semilogy(10:10:50,abs(rek-corr),'r-*',10:10:50,abs(expl-corr),'g-d')
title('computation error')
xlabel('n')
ylabel('error')
legend('recursive','explicit',-1)

print -dps 'ha1.ps'

How can I do that?

Answer 1

You are on the right track. The recursion can be written as

function I=intrek(n,I0)
  I = I0;
  for i = 1:n
    I = e-i*I;
  endfor
endfunction

Or, if you prefer a while loop

function I=intrek(n,I0)
  I = I0;
  i = 0;
  while i < n
    i++;
    I = e-i*I;
  endwhile
endfunction

The summation can be written as

function I = intexp(n)
  I = factorial(n)*(e*sum((-1).^(0:n) ./ factorial(n - (0:n))) - (-1)^n);
endfunction

You can also find the integral using general-purpose numeric integration routines as

In = @(n) integral(@(x) x .^ n .* exp(x), 0, 1);