# Time-stepping an equation in Mathematica

I have an iterative integral equation of the form $$n(k,t+dt)=n(k,0)+dt\int_{0}^{\infty}dp\int_{0}^{\pi}d\alpha\int_{0}^{2\pi}d\beta\int_{0}^{\pi}da\,p^{2}q^{2}\sin(\alpha)\sin(a)\frac{(n(p,t)+1)}{\omega(p)}$$ where $\omega(p) =p^{2}+1$, and the value of $\lvert q\rvert=Q$ is set by the equation $\sqrt{(k-p-q)^{2}}=0$. There is also the initial condition that $n(p,0)=0$.

I want to write a code that constructs an array of values of $n(k,t)$ for each value of momentum $k$ at each instant in time $t$. I need the code to time-step the above equation, such that for each subsequent time-step, $t+dt$ the code uses the value of $n$ from the previous time-step, i.e. $n(p,t)$, in order to determine the value of $n$ at $t+dt$, i.e. $n(k,t+dt)$.

My strategy so far has been to discretise the integral over $p$ and calculate as follows:

p = {P*Sin[\[Alpha]]*Cos[\[Beta]], P*Sin[\  [Alpha]]*Sin[\[Beta]],P*Cos[\[Alpha]]};
k = {0, 0, X};
q = {Q*Sin[a], 0, Q*Cos[a]};
q = q/.Q->Solve[Sqrt[(k - p - q).(k - p - q) + 1] == 0, Q][[1,1]];
omega[p_, q_, k_] := Sqrt[(k - p - q).(k - p - q) + 1]
dP=0.1;
n[X_, m*dt] = #1 + m*dt*2*Pi*Sum[dP*Integrate[(p.p)*(q.q)*(1/omega[p,0,k])*Sin[\[Alpha]]*Sin[a]*(#2 + 1), {\[Alpha], 0, Pi}, {\[Beta], 0, 2*Pi}, {a, 0, Pi}], {P, 0, 100, 0.1}] &
dt = 0.1;
n[k_, 0] = 0.0;
Array[n[X, m*dt] /@ {n[k, 0], n[k-p-q, (m - 1)*dt]}, {X, 0, 100}, {m, 0, 100}]


such that, starting at $t=0$, the array contains values for $n$, incremented by an amount $dt$ in each subsequent row (up to $n(X,mdt)$). However, what I written above doesn't seem quite right as it stands. [I'm fairly new to the concept of pure functions and maps in Mathematica.]

Any suggestions of how to proceed would be much appreciated.

• Two things: 1) How does $F$ change from $t$ to $t+dt$? You only show how to update $n$. 2) With the initial condition you have posted, it looks like the integrand is 0 every time, so you will get $n=0$ for all $t$. – Marius Ladegård Meyer Nov 10 '16 at 17:30
• @MariusLadegårdMeyer I've updated the question to include a less general function, that will suffice as an example of how to proceed (the actual case I need to look at is a bit more complicated, but not too far from this). – user35305 Nov 10 '16 at 17:36
• Hm, can you explain how the Mathematica code you have added relates to the first equation...? There is no q, omega etc. in the equation. And again, with $n(k, 0) = 0$, how do you expect to get anywhere? Choosing $t=0$ in the equation gives $n(k, dt) = 0$ and the same for all $m dt$, $m = 1,2,...$. – Marius Ladegård Meyer Nov 10 '16 at 21:48
• @MariusLadegårdMeyer Apologies, I was thinking it through in my head and forgot to write the extra details in the end (it was very sloppy of me, sorry). Have added them now. In principle, it should be solvable right, since at $n(p,0)$ the integrand is non-zero?! – user35305 Nov 10 '16 at 22:55
• Sorry, I still have some uncertainties that I need to clarify in order to be able to write a good answer...In the text you say $q$ is found by solving $(k - p - q)^2 = 0$, but in the code you solve Sqrt[(k - p - q).(k - p - q) + 1] == 0 which is not the same. Then you set omega to be the value of Sqrt[(k - p - q).(k - p - q) + 1] at the q solved for above, but that solution MUST give omega equal to 0! And you cannot divide by 0 in the integrand.. And later you set x equal to (k - p - q).(k - p - q) at the q value solved for, but that MUST give -1, right? – Marius Ladegård Meyer Nov 12 '16 at 18:11