1
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.

Michael E2
  • 235,386
  • 17
  • 334
  • 747
Rupesh
  • 887
  • 5
  • 10
  • 1
    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)? – Szabolcs Mar 11 '19 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. – Henrik Schumacher Mar 11 '19 at 14:35
  • Try e.g. Rx[1, 2, t], and notice that the integral is not solved. – Szabolcs Mar 11 '19 at 14:38
  • 1
    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). – Szabolcs Mar 11 '19 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? – Rupesh Mar 11 '19 at 14:45
  • thanks henrik I will go through that too – Rupesh Mar 11 '19 at 14:49
  • Use NIntegrate like this: Rx[x_?NumericQ, y_?NumericQ, t_?NumericQ] := Dx[x, y, t] E^(-I NIntegrate[AA[x, y, r], {r, -3, t}]) – Michael E2 Mar 11 '19 at 17:32
  • Using a slightly higher order integration rule saves a couple of seconds: NIntegrate[AA[x, y, r], {r, -3, t}, Method -> {"GaussKronrodRule", "Points" -> 11}] – Michael E2 Mar 11 '19 at 17:34
  • thanks michael i will try it – Rupesh Mar 11 '19 at 17:36

2 Answers2

2

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?

Michael E2
  • 235,386
  • 17
  • 334
  • 747
  • Thanks Michael, I will give it a go the way you said. – Rupesh Mar 11 '19 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. – Michael E2 Mar 11 '19 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}] – Rupesh Mar 11 '19 at 19:18
1

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.

Michael E2
  • 235,386
  • 17
  • 334
  • 747
  • It seems to work. I surely need to work on the details of the code. Thank you for your time and help. – Rupesh Mar 12 '19 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. – Rupesh Mar 13 '19 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? – Michael E2 Mar 13 '19 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.) – Rupesh Mar 13 '19 at 20:51
  • 1
    @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. – Michael E2 Mar 14 '19 at 00:34
  • Thanks Michael, you have been a great help. It makes sense. – Rupesh Mar 14 '19 at 14:09