# Slow processing of coupled differential equation using NDsolve

ClearAll["Global*"]
AA[x_, y_, t_] := Sqrt[3 + 2*Cos[(t*x)/E^(t^2/2)] + Cos[2*y] + 2*Cos[(t*y)/E^(t^2/2)]]

Dx[x_, y_,t_] := (Cos[2*(3/E^9 + t/E^t^2 + x)] + 2*Cos[3*(3/E^9 + t/E^t^2) + x]*Cos[y] + (2*I)*Cos[y]*Sin[2*(x + (3/E^9 + t/E^t^2)*x)])/ Sqrt[1 + Abs[Cos[2*(3*(3/E^9 + t/E^t^2) + x)]]]
F[t_] :=  E^-t^2 (1 - 2 t^2)
Rx[x_, y_, t_] := Dx[x, y, t] E^(-I Integrate[AA[x, y, r], {r, -3, t}])

L[x_, y_, t_] := -I F[t] Rx[x, y, t]
M[x_, y_, t_] := -I F[t] Conjugate[Rx[x, y, t]]

eqn[x_, y_, t_] := {A1'[t] == L[x, y, t] A2[t],
A2'[t] == A1[t] M[x, y, t],
A1[-3] == 0, A2[-3] == 1}
sol := {A1, A2} /. (NDSolve[eqn[x, y, t],{A1, A2}, {t, -3, 3}])
P=Table[sol, {x, 1, 2}, {y, 1, 2}]


I am having a problem getting a solution. It just keeps running. And, considering the fact I need it for a wide range of x and y, I definitely need some help.

• Try to construct a minimal example. Are you sure it's NDSolve that takes long? I think it's Integrate. Have you tried running parts of your code individually (which should always be the first debugging step you take)? Mar 11, 2019 at 14:35
• NDSolve has a hard time to solve a differential equation with symbolic parameters x and y. Have a look into ParametricNDSolve which, of course, includes reading its documentation. Mar 11, 2019 at 14:35
• Try e.g. Rx[1, 2, t], and notice that the integral is not solved. Mar 11, 2019 at 14:38
• I suggest rephrasing the problem so that it does not include an integral. This is often possible by introducing an additional function to solve for in NDSolve, representing the integral (differentiating it by t just gives the integrand). Mar 11, 2019 at 14:41
• thanks szabolcs, as you said the problem comes from the integral in Rx. Is there any way I can get past it? Mar 11, 2019 at 14:45

Using NDSolve to compute the indefinite integral saves a lot time, when the integral will be evaluated many times as in this appliction:

ClearAll[int];
mem : int[x_, y_] := (* mem represents the function call and "memoizes" the result *)
mem = Block[{i, r, t},
NDSolveValue[{i'[r] == AA[x, y, r], i[-3] == 0}, i, {r, -3, 3},
PrecisionGoal -> 12,  (* integrand is nice: squeeze out a bit more precision *)
InterpolationOrder -> All]];  (* ensure accuracy between interpolation nodes *)
Rx[x_?NumericQ, y_?NumericQ, t_?NumericQ] :=
Dx[x, y, t] E^(-I int[x, y][t]);


This brings the execution time down to around 0.15 sec. (Simply replacing Integrate with NIntegrate brought the time down to around 7 sec.)

Here is an explanation of memoization: What does the construct f[x_] := f[x] = ... mean?

• Thanks Michael, I will give it a go the way you said. Mar 11, 2019 at 19:03
• @maeinss You're welcome. I think Szabolcs' idea is a bit better, but I don't know if he will have time to answer. Mar 11, 2019 at 19:08
• that was amazing I never thought it that way. However, even with the value of Rx the solution still keeps running. I am taking about last few steps. eqn[x_, y_, t_] := {A1'[t] == L[x, y, t] A2[t],A2'[t] == A1[t] M[x, y, t],A1[-3] == 0, A2[-3] == 1} sol := {A1, A2} /. (NDSolve[eqn[x, y, t],{A1, A2}, {t, -3, 3}]) P=Table[sol, {x, 1, 2}, {y, 1, 2}] Mar 11, 2019 at 19:18

Here is one interpretation of Szabolcs's idea, using the variable i[t] for the integral:

ClearAll[Rx, eqn];
Rx[x_, y_, t_] := Dx[x, y, t] E^(-I i[t]);  (* i[t] is the integral up to t *)

eqn[x_, y_, t_] = {A1'[t] == L[x, y, t] A2[t],
A2'[t] == A1[t] M[x, y, t], A1[-3] == 0, A2[-3] == 1,
i'[t] == AA[x, y, t], i[-3] == 0};  (* the last two eqns define the integral *)

PrintTemporary@Dynamic@{x, y, Clock[Infinity]};
sol := {A1, A2} /. (NDSolve[eqn[x, y, t], {A1, A2}, {t, -3, 3}])
P = Table[sol, {x, 1, 2}, {y, 1, 2}] // AbsoluteTiming


It takes around 0.03 sec., significantly faster than my (other) answer.

• It seems to work. I surely need to work on the details of the code. Thank you for your time and help. Mar 12, 2019 at 14:06
• Hi Michael, I was analyzing the code and I have one confusion. How do you specify i[-3]=0? Integral of AA[x, y, t] at t=-3 can have any value. If so, how can we evaluate. Mar 13, 2019 at 18:05
• @maeinss I specified i[-3] with the code i[-3] == 0; that's how I did it. Perhaps you're asking why? i[t] represents Integrate[AA[x, y, r], {r, -3, t}]. For t == -3, the integral is zero, no matter what the value of AA[x, y, t], since the end points of integration are both r == -3.. Or are you asking something else entirely? Mar 13, 2019 at 18:32
• Sorry, I stated it wrong, I meant why. I can't see a reason why integration of AA[x,y,r/t] should be 0 at the boundary. A1[-3] == 0, A2[-3] == 1 because they are probabilities. However, even at those endpoint values system can have finite AA value other than 0 and so its integration can be non-zero. (In my case AA is the Energy.) Mar 13, 2019 at 20:51
• @maeinss It's basic calculus: $\int_{-3}^{-3} f(r)\;dr = 0$ by definition, no matter what $f(r)$ is (provided $f$ is defined at $-3$). Whether $f(-3) = 10$ or $f(-3) = -100$, the integral is zero. Compare by executing Integrate[AA[x, y, r], {r, -3, -3}]` in Mathematica. Mar 14, 2019 at 0:34