Mathematica 8.0  solve two recursive relationsby EigenCake Tags: mathematica, recursive, relations, solve 

#1
Nov1212, 01:01 PM

P: 13

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 buildin function in Mathematica), however, I am very new in using Mathematica. Please help! Thank you~ 



#2
Nov1312, 12:31 PM

P: 972

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[n1,x] Sum[e[j]y[nj,x],{j,1,n1}]), {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](w[s]y[n1,s]Sum[e[j]y[nj,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/2x^2/2)(3+x^4) Erfi[x/Sqrt[2]] does not converge on {∞, ∞}. From In[7]:= Integrate:: idiv: Integral of E^(t^2/2x^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)} 



#3
Nov1412, 02:25 AM

P: 13

Wow, Bill, you are amazing!
So far as I know, there is a mistake on your code On you code:
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:
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 wanna know whether your code is able to give right answer or not. Thank you so much, Bill! 



#4
Nov1412, 02:48 AM

P: 13

Mathematica 8.0  solve two recursive relations
Now I know how to run your code on my computer
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 halffreezed. 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! :) 



#5
Nov1412, 01:11 PM

P: 972

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[n1,x]Sum[e[j] y[nj,x],{j,1,n1}]),{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[n1,s]Sum[e[j] y[nj,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[n1,x]Sum[e[j] y[nj,x],{j,1,n1}]),{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[s] y[n1,s]Sum[e[j] y[nj,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. 



#6
Nov1412, 10:58 PM

P: 13

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! 


Register to reply 
Related Discussions  
Mathematica: Need help to solve this pde set  Math & Science Software  2  
finding recursive relations in Mathematica  Math & Science Software  4  
how do you solve third order recurrence relations?  Calculus & Beyond Homework  1  
Recursive Function in Mathematica  Math & Science Software  10  
Recursive Relations.  Set Theory, Logic, Probability, Statistics  0 