First, you might recognize that 1/e = 1/exp(1) = exp(-1).
So if you are going to compute an exponential, thus exp(1), then why not save time and just compute exp(-1) in the first place?????
Next, do NOT call a function integral!!!!!!!!!!!!! If you do, then the existing, VERY useful function named integral will no longer be accessible. And then your next frenzied question will be about why integral does not work anymore. So call your function something like myintegral instead.
Next, your real question is what do you do with y(n-1). But you already have a formula for y(n). This is a basic recurrence relation. Just call the function recursively!
By the way, you don't need a return statement in the if statement, if your code will just exit the function anyway.
Finally, you set the result to y(n). But your function takes a value for n, and then returns a SCALAR value for y_n. However, if you try to stuff it into y(n), then you will crete the VECTOR y, that is n elements long. All of the other elements will be zero, except the nth element of that vector.
function y = myint1(n)
%INTEGRAL
if n == 1
y = 1 - 2*exp(-1);
end
y = n .* myint1(n-1) - exp(-1);
end
Of course, this is slightly inefficient if n would be at all large, because it will require you to call the function integral n times.
Or, perhaps, what you really want is to return the value for all values of y(n) up to but not exceeding n. Then it would be efficient to just write it as a loop.
function y = myint2(n)
%INTEGRAL
y = zeros(1,n);
expminus1 = exp(-1);
y(1) = 1 - 2*expminus1;
for m = 2:n
y(m) = m.*y(m-1) -expminus1;
end
That code would return a vector y, for all values of up to n.
With n==20, lets see how the latter code works.
format long g
>> y
y =
Columns 1 through 5
0.264241117657115 0.160602794142788 0.113928941256923 0.0878363238562483 0.0713021781097991
Columns 6 through 10
0.0599336274873523 0.0516559512400239 0.0453681687487486 0.0404340775672951 0.0364613345015088
Columns 11 through 15
0.0331952383451544 0.03046341897041 0.0281450054438883 0.0261506350429941 0.024380084473469
Columns 16 through 20
0.0222019104040609 0.00955303569759369 -0.195924798614756 -4.0904506148518 -82.1768917382075
What is interesting is that it starts to yield negative numbers for n as large as 20. Even 18 seems to be a problem. (By the way, note that I actually need to use the true functino integral here to compare results!)
m = 1;[m,y(m),integral(@(x) kernel(x,m),0,1)]
ans =
1 0.264241117657115 0.264241117657115
>> m = m+1;[m,y(m),integral(@(x) kernel(x,m),0,1)]
ans =
2 0.160602794142788 0.160602794142788
>> m = m+1;[m,y(m),integral(@(x) kernel(x,m),0,1)]
ans =
3 0.113928941256923 0.113928941256923
>> m = m+1;[m,y(m),integral(@(x) kernel(x,m),0,1)]
ans =
4 0.0878363238562483 0.0878363238562491
>> m = m+1;[m,y(m),integral(@(x) kernel(x,m),0,1)]
ans =
5 0.0713021781097991 0.0713021781098032
>> m = m+1;[m,y(m),integral(@(x) kernel(x,m),0,1)]
ans =
6 0.0599336274873523 0.0599336274873766
>> m = m+1;[m,y(m),integral(@(x) kernel(x,m),0,1)]
ans =
7 0.0516559512400239 0.0516559512401941
>> m = m+1;[m,y(m),integral(@(x) kernel(x,m),0,1)]
ans =
8 0.0453681687487486 0.0453681687501108
>> m = m+1;[m,y(m),integral(@(x) kernel(x,m),0,1)]
ans =
9 0.0404340775672951 0.040434077579555
>> m = m+1;[m,y(m),integral(@(x) kernel(x,m),0,1)]
ans =
10 0.0364613345015088 0.0364613346241073
>> m = m+1;[m,y(m),integral(@(x) kernel(x,m),0,1)]
ans =
11 0.0331952383451544 0.0331952396937377
>> m = m+1;[m,y(m),integral(@(x) kernel(x,m),0,1)]
ans =
12 0.03046341897041 0.0304634351534098
>> m = m+1;[m,y(m),integral(@(x) kernel(x,m),0,1)]
ans =
13 0.0281450054438883 0.0281452158228847
>> m = m+1;[m,y(m),integral(@(x) kernel(x,m),0,1)]
ans =
14 0.0261506350429941 0.026153580348944
>> m = m+1;[m,y(m),integral(@(x) kernel(x,m),0,1)]
ans =
15 0.024380084473469 0.0244242640627179
>> m = m+1;[m,y(m),integral(@(x) kernel(x,m),0,1)]
ans =
16 0.0222019104040609 0.0229087838320443
>> m = m+1;[m,y(m),integral(@(x) kernel(x,m),0,1)]
ans =
17 0.00955303569759369 0.0215698839733108
>> m = m+1;[m,y(m),integral(@(x) kernel(x,m),0,1)]
ans =
18 -0.195924798614756 0.0203784703481514
Look carefully at what you see here. The last few digits start to go bad as n gets larger, and it does not even require you to go very far. Even at m==4, we see some deviation. And one thing that we should recognize is that the integral in question should NEVER yield a negative number, since all terms are positive, over the interval [0,1]. So how could the integral be negative?
What is happening? Why does this fail?
You need to see that for larger values of n, we are multiplying by n in this recursive relation. And we are also subtracting two numbers.
n .* myint1(n-1) - exp(-1);
So what happens is any tiny error in y(n-1) gets amplified, because you just multiplied it by n. Think of it like the factorial function. nd what is factorial(18), or a number that large?
factorial(18)
ans =
6.402373705728e+15
So, if there is an error in the least significant bit of y(1), by only one bit in that decimal place, then it gets multiplied (AMPLIFIED!) by something on the order of 6e15 by the time you get to y(18).
Hmm. Can we fix this? Well, yes. There is a trick we can use, and it is one that is often of value in problems like this. The ide is to start out past the END, and work backwards.
That is, if we wish to compute the value of y(n), start at n+20. We have this recursion:
y(n) = n .* y(n-1) - exp(-1);
Rewrite it as:
y(n-1) = (y(n) + exp(-1))/n
Now we can write it in a loop, going backwards.
function yn = myint3(n)
%INTEGRAL, starting past the top...
yn = .1; % Anything even close will suffice here.
N = n + 20;
expminus1 = exp(-1);
for m = N:-1:(n+1)
yn = (yn + expminus1)/m;
end
Does this work?
m=1,integral(@(x) kernel(x,m),0,1),myint3(m)
m =
1
ans =
0.264241117657115
ans =
0.264241117657115
>> m=18,integral(@(x) kernel(x,m),0,1),myint3(m)
m =
18
ans =
0.0203784703481514
ans =
0.0203784703481514
>> m=30,integral(@(x) kernel(x,m),0,1),myint3(m)
m =
30
ans =
0.0122495029578589
ans =
0.0122495029578589
As you see, now this code with the backwards loop works for virtually any value of n, whereas when you use the forward loop, it starts to fail very early. We can even compare it to a symbolic integral.
m = 18
m =
18
>> int(X.^m*exp(-X),0,1)
ans =
6402373705728000 - 17403456103284421*exp(-1)
>> vpa(ans)
ans =
0.020378470348151434104822774864287
>> myint3(18)
ans =
0.0203784703481514
Again, it yields essentially full precision results. Why did it work? Think about the difference in terms of how I formulated the recurrence. Dividing by large numbers will decrease any noise, until it all gets killed off, even if I started with complete crap for my initial guess at y(n+20).
This is a good trick to remember. (You may forget it by tomorrow, but perhaps not, or perhaps it will ring a bell in some later class.)
