Mathematica Mathematica 8.0 - solve two recursive relations

[EigenCake]

Hi, all!

I have trouble by using Mathemtica to solve the following problem as shown on the attachment. You see that the two recursive relations depend one another. I plan to write a "For" loop to evaluate E, instead of using Rsolve (a build-in function in Mathematica), however, I am very new in using Mathematica. Please help! Thank you~

 

[Bill Simpson]

Check this carefully, perhaps there is an error that can be corrected.

In[1]:= e[0]=1/2;
y[0,x_]:=E^-(x^2/4);
w[x_]:=x^4/4;
e[n_]:=Integrate[y[0,x](w[x]y[n-1,x]- Sum[e[j]y[n-j,x],{j,1,n-1}]), {x,-Infinity,Infinity}]/ Integrate[y[0,x]^2,{x,-Infinity,Infinity}];
y[n_,x_]:=y[0,x]Integrate[1/y[0,t]^2,{t,0,x}] Integrate[y[0,s](wy[n-1,s]-Sum[e[j]y[n-j,s],{j,1,n}]), {s,-Infinity,t}];

In[6]:= {e[1],y[1,x]}

Out[6]= {3/4, -(E^(-t^2/2 - x^2/4)*Sqrt[Pi/2]*t*(3 + t^2)*Erfi[x/Sqrt[2]])/4}

In[7]:= {e[2],y[2,x]}

From In[7]:= Integrate:: idiv: Integral of E^-(t^2/2-x^2/2)(-3+x^4) Erfi[x/Sqrt[2]] does not converge on {-∞, ∞}.

From In[7]:= Integrate:: idiv: Integral of E^-(t^2/2-x^2/2)(-3+x^4) Erfi[x/Sqrt[2]] does not converge on {-∞, ∞}.

From In[7]:=Power:: infy: Infinite expression 1/0 encountered.

From In[7]:=Power:: infy: Infinite expression 1/0 encountered.

From In[7]:=Power:: infy: Infinite expression 1/0 encountered.

From In[7]:= General:: stop: Further output of Power :: infy will be suppressed during this calculation.

Out[7]= {Integrate[((3*E^(-t^2/2 - x^2/4)*Sqrt[Pi/2]*t*(3 + t^2)*Erfi[x/Sqrt[2]])/16 - (E^(-t^2/2 - x^2/4)*Sqrt[Pi/2]*t*(3 + t^2)*x^4*Erfi[x/Sqrt[2]])/16)/E^(x^2/4), {x, -Infinity, Infinity}]/Sqrt[2*Pi], (Sqrt[Pi/2]*Erfi[x/Sqrt[2]]*Integrate[((3*E^(-s^2/4 - t^2/2)*Sqrt[Pi/2]*t*(3 + t^2)*Erfi[s/Sqrt[2]])/16 - (E^(-s^2/4 - t^2/2)*Sqrt[Pi/2]*s^4*t*(3 + t^2)*Erfi[s/Sqrt[2]])/16 - Integrate[((3*E^(-t^2/2 - x^2/4)*Sqrt[Pi/2]*t*(3 + t^2)*Erfi[x/Sqrt[2]])/16 - (E^(-t^2/2 - x^2/4)*Sqrt[Pi/2]*t*
(3 + t^2)*x^4*Erfi[x/Sqrt[2]])/16)/E^(x^2/4), {x, -Infinity, Infinity}]/(E^(s^2/4)*Sqrt[2*Pi]))/E^(s^2/4), {s, -Infinity, t}])/E^(x^2/4)}

 

[EigenCake]

Wow, Bill, you are amazing!

So far as I know, there is a mistake on your code

On you code:

y[n_,x_]:=y[0,x]Integrate[1/y[0,t]^2,{t,0,x}] Integrate[y[0,s](w[s]y[n-1,s]-Sum[e[j]y[n-j,s],{j,1,n}]), {s,-Infinity,t}];

Because the upper bound of the integral "ds" is "t", so the integral "dt" includes the integration over "ds", which means here is to evaluate a double integral, instead of evaluating two separate and independent integrals. So I think it must be:

y[n_,x_]:=y[0,x]Integrate[(1/y[0,t]^2)*Integrate[y[0,s](w[s]y[n-1,s]-Sum[e[j]y[n-j,s],{j,1,n}]), {s,-Infinity,t}], {t, 0, x}];

I don't know whether the modified code I gave to you will give correct answer or not, because until now I do NOT know how to run your code on my computer properly. Whatever, if the above code does not give the correct answer, which is:
y1 = (-((3*x^2)/8) - x^4/16)/E^(x^2/4)
then the code is wrong.

Before I post my message, I wrote the following code to calculate y1, by properly setting value for y0[t], then the following code:

y1[x_] == y0[x]*Integrate[(1/((y0[t])^2))*Integrate[y0[s]*((s^4/4)*y0[s] - (3/4)*y0[s]), {s, -Infinity, t}], {t, 0, x}]

gives Output:

(1/8 E^(-(x^2/
     4)) (-3 + x^4) (\[Pi] Erfi[x/Sqrt[2]] + 
      x^2 HypergeometricPFQ[{1, 1}, {3/2, 2}, x^2/2]))[x_] == 
 E^(-(x^2/4)) (-((3 x^2)/8) - x^4/16)

By observing this Output, the last part E^(-(x^2/4)) (-((3 x^2)/8) - x^4/16) is the right answer, but I don't know why Mathematica also give an expression about HypergeometricPFQ, which is not only useless but also ruins the result to calculating E2, so that E2 is just a piece of mess!

If use y1 = (-((3*x^2)/8) - x^4/16)/E^(x^2/4) to calculate E2, then E2 = -21/8.

The correct answer for the first four E's are:
E0 = 1/2;
E1 = 3/4;
E2 = -21/8;
E3 = 333/16;

If any code does give output E's like those, then something must be wrong in the code.

Okay, now I need to figure out a way to run your code on my computer. Could you correct that mistake on your code please? After correction, I want to know whether your code is able to give right answer or not.

Thank you so much, Bill!

 

[EigenCake]

Now I know how to run your code on my computer

in[1] := e[0] = 1/2;
y[0, x_] := Exp[-(x^2/4)];
w[x_] := x^4/4;
e[n_] := Integrate[
    y[0, x] (w[x] y[n - 1, x] - 
       Evaluate[Sum[e[j] y[n - j, x], {j, 1, n - 1}]]), {x, -Infinity,
      Infinity}]/Integrate[y[0, x]^2, {x, -Infinity, Infinity}];
y[n_, x_] := 
  y[0, x] Integrate[(1/y[0, t]^2)*
     Integrate[
      y[0, s] (w[s] y[n - 1, s] - 
         Evaluate[Sum[e[j] y[n - j, s], {j, 1, n}]]), {s, -Infinity, 
       t}], {t, 0, x}] ;

This code gives correct answers for all four E's mentioned on my last message. So I believe that the code is able to give correct answer to all E's.

However, the computation is so slow that I have to wait for 10 minutes in order to have the answer E5 = 916731/256, and my computer is half-freezed. This is slower than what I expected. So I have to find another way to find E200, instead of using the recursive relations.

Nonetheless, the problem on the main message is solved, and the code is quite neat! Thank you a lot, Bill! I highly appreciate your help! You're amazing! :)

 

[Bill Simpson]

Now that you have it (hopefully) correct, let's speed it up.

First your latest modified code, unchanged, hopefully this is correct. I did remove the Evaluate[] when I realized earlier those were not necessary.

In[1]:= e[0]=1/2;
y[0,x_]:=Exp[-(x^2/4)];
w[x_]:=x^4/4;
e[n_]:=Integrate[y[0,x] (w[x] y[n-1,x]-Sum[e[j] y[n-j,x],{j,1,n-1}]),{x,-Infinity,Infinity}]/Integrate[y[0,x]^2,{x,-Infinity,Infinity}];
y[n_,x_]:=y[0,x] Integrate[(1/y[0,t]^2)*Integrate[y[0,s] (w y[n-1,s]-Sum[e[j] y[n-j,s],{j,1,n}]),{s,-Infinity,t}],{t,0,x}];

In[6]:= Timing[{y[5,x],e[5]}]

Out[6]= {450.457*Second, {((-916731*x^2)/512 - (651363*x^4)/2048 - (89673*x^6)/2048 - (69107*x^8)/16384 - (250183*x^10)/786432 - (29177*x^12)/1572864 - (1325*x^14)/1572864 - (269*x^16)/9437184 - (25*x^18)/37748736 - x^20/125829120)/E^(x^2/4), 916731/256}}

Now kill the kernel and restart with a tiny change to your code

In[1]:= e[0]=1/2;
y[0,x_]:=Exp[-(x^2/4)];
w[x_]:=x^4/4;
e[n_]:=e[n]=Integrate[y[0,x] (w[x] y[n-1,x]-Sum[e[j] y[n-j,x],{j,1,n-1}]),{x,-Infinity,Infinity}]/Integrate[y[0,x]^2,{x,-Infinity,Infinity}];
y[n_,x_]:=y[n,x]=y[0,x] Integrate[(1/y[0,t]^2)*Integrate[y[0,s] (w y[n-1,s]-Sum[e[j] y[n-j,s],{j,1,n}]),{s,-Infinity,t}],{t,0,x}];

In[6]:= Timing[{y[5,x],e[5]}]

Out[6]= {12.628*Second, {((-916731*x^2)/512 - (651363*x^4)/2048 - (89673*x^6)/2048 - (69107*x^8)/16384 - (250183*x^10)/786432 - (29177*x^12)/1572864 - (1325*x^14)/1572864 - (269*x^16)/9437184 - (25*x^18)/37748736 - x^20/125829120)/E^(x^2/4), 916731/256}}

That is a very tiny change that makes this 35 times faster, and will make it far far far far faster if you are going to try n=200, that should not break the mathematics. It is just a Mathematica trick. I'm not even going to try to explain it.

As always, check this, test this, verify this is correct. It is way too easy to type a bunch of stuff in, get something out that isn't obviously completely wrong and think you are done.

Note: If you are going to use Mathematica, particularly if you are going to attack big problems then you should seriously consider getting the fastest CPU and memory that you can afford. A good recent fairly fast I5 or I7 or a comparable AMD CPU along with 16 or 32 gig of memory, all dedicated to the sole task of running Mathematica, will make your life much easier. But, even with that, some problems are hard and take a long time, even if you have written the best code implementing the best algorithm that you can, and it is not uncommon to wait hours or even in extreme cases weeks or months to get a solution for a really hard problem.

 

[EigenCake]

Oh my god! You save me lots of effort!

I tried your code and get E5 and y5 in 10.407 second. That is so awesome! I don't know such Mathematica trick, and the reason behind it and how it works, but soon or later I will know it, after manipulating Mathematica for a while.

Now this problem seems COMPLETELY solved, and you are the MVP! Thank you so much, Bill~~~ And have a wonderful day!