1
$\begingroup$ 
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.

Improve this question
$\endgroup$
9
  • 1
    $\begingroup$ 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)? $\endgroup$
    – Szabolcs
    Mar 11, 2019 at 14:35
  • $\begingroup$ 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. $\endgroup$
    – Henrik Schumacher
    Mar 11, 2019 at 14:35
  • $\begingroup$ Try e.g. Rx[1, 2, t], and notice that the integral is not solved. $\endgroup$
    – Szabolcs
    Mar 11, 2019 at 14:38
  • 1
    $\begingroup$ 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). $\endgroup$
    – Szabolcs
    Mar 11, 2019 at 14:41
  • $\begingroup$ thanks szabolcs, as you said the problem comes from the integral in Rx. Is there any way I can get past it? $\endgroup$
    – Rupesh
    Mar 11, 2019 at 14:45

2 Answers 2

Reset to default
2
$\begingroup$

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?

Improve this answer
$\endgroup$
3
  • $\begingroup$ Thanks Michael, I will give it a go the way you said. $\endgroup$
    – Rupesh
    Mar 11, 2019 at 19:03
  • $\begingroup$ @maeinss You're welcome. I think Szabolcs' idea is a bit better, but I don't know if he will have time to answer. $\endgroup$
    – Michael E2
    Mar 11, 2019 at 19:08
  • $\begingroup$ 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}] $\endgroup$
    – Rupesh
    Mar 11, 2019 at 19:18
1
$\begingroup$

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.

Improve this answer
$\endgroup$
6
  • $\begingroup$ It seems to work. I surely need to work on the details of the code. Thank you for your time and help. $\endgroup$
    – Rupesh
    Mar 12, 2019 at 14:06
  • $\begingroup$ 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. $\endgroup$
    – Rupesh
    Mar 13, 2019 at 18:05
  • $\begingroup$ @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? $\endgroup$
    – Michael E2
    Mar 13, 2019 at 18:32
  • $\begingroup$ 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.) $\endgroup$
    – Rupesh
    Mar 13, 2019 at 20:51
  • 1
    $\begingroup$ @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. $\endgroup$
    – Michael E2
    Mar 14, 2019 at 0:34

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.