1
$\begingroup$

I am trying to plot results for a very simple differential equation of the form:

$$\frac{\partial^2 x(N,z'(N))}{\partial N^2} = F(N,z'(N)), $$

where $z'(N)$ is a function of $N$ that needs to be solved using FindRoot for at every $N$ position, and $F(N,z')$ is a nasty equation that results from a numerical integration over:

$$ F(N,z') = \int_{-\infty}^{\infty} \exp\left( -\frac{x'^2}{2\sigma_{x'}^2} \right) F(N,z',x')dx'$$.

So, I put together some mathematica code but it runs terrible slow (on the order of a day or two)! I noticed there were some things that affected the speed of the code particularly the numerical coefficient in front of $F(N,z'(N))$. But I was wondering if there is any help to be given to get better/faster results! Any help would be greatly appreciated!

Note: I had to use $NN$ in place of $N$ in my equations because its a function in mathematica. Also, in the FN function I have to actually copy and paste the output of FNzprime (an ugly mess) in the integrand for it to evaluated.

(*constants*)
e = -1.60217733*10^-19;
m = 9.109389699999999*10^-31;
epsilon = 8.854187817620391*10^-12;

(*basic equations*)
rs2 = {zprime, xprime + K/(gamma*kw) Sin[kw*zprime], 0};
ro2 = {(NN + 10000)*lw, x + K/(gamma*kw) Sin[kw*(NN + 10000)*lw], 0};

betas = {beta - K^2/(4 gamma^2) Cos[2 kw*zprime],K/gamma Cos[kw*zprime], 0};
betao = {beta - K^2/(4 gamma^2) Cos[2 kw*(NN + 10000)*lw],K/gamma Cos[kw*(NN + 10000)*lw], 0};

betaDot = {(c*K^2*kw)/(2 gamma^2)Sin[2 kw*zprime], -((c*K*kw)/gamma) Sin[kw*zprime], 0};

deltar2 = ro2 - rs2;
Rgam2 = Sqrt[deltar2[[1]]^2 + deltar2[[2]]^2];
Ec2 = (e/(4 Pi*epsilon)) (deltar2/Rgam2 - betas)/(gamma^2 Rgam2^2 (1 - (deltar2/Rgam2).betas)^3);
Erad2 = (e/(4 Pi*epsilon)) Cross[deltar2/Rgam2, Cross[deltar2/Rgam2 - betas, betaDot]]/(c*Rgam2*(1 - (deltar2/Rgam2).betas)^3);

Bc2 = Cross[deltar2/Rgam2, Ec2];
Brad2 = Cross[deltar2/Rgam2, Erad2];

Fbc2 = Cross[betao, Bc2];
Fbrad2 = Cross[betao, Brad2];


sumEtran = (Ec2[[2]] + Erad2[[2]]);
sumFBtran = Fbc2[[2]] + Fbrad2[[2]];



(*Numeric Functions*)

ZPRIME[NN_?NumericQ, x_?NumericQ, xprime_?NumericQ, gamma_, K_, kw_, beta_, sigma_, lw_] :=zprime /. FindRoot[sigma == (1/(gamma kw))Sqrt[gamma^2 + K^2] (EllipticE[kw*(NN + 10000)*lw, K^2/(gamma^2 + K^2)] - EllipticE[kw zprime, K^2/(gamma^2 + K^2)]) - beta \[Sqrt](((NN + 10000)*lw - zprime)^2 + (x - xprime + (K Sin[kw *(NN + 10000)*lw])/(gamma kw) - (K Sin[kw zprime])/(gamma kw))^2), {zprime, 0}]


coeff = ((e*lw^2)/(gamma*m*beta^2*c^2) (10^-10/e)/(2 Pi (30*10^-6) (10^-5)) Exp[-(sigma^2/(2 (10^-5)^2))]);
FNzprime =coeff (sumEtran + sumFBtran) /. {lw -> 0.026, K -> 1, beta -> Sqrt[1 - 1/(4000/0.511)^2], gamma -> 4000/0.511, c -> 3*10^8, kw -> (2 Pi)/0.026, zprime -> ZPRIME}

FN[NN_?NumericQ, x_?NumericQ, sigma_?NumericQ] :=With[{ZPRIME = ZPRIME[NN, x, 0, 4000/0.511, 1, (2 Pi)/0.026, Sqrt[1 - 1/(4000/0.511)^2], sigma, 0.026]}, 
  NIntegrate[ (Exp[-(xprime^2/(2 (30*10^-6)^2))]) FNzprime, {xprime, -300*10^-6, 300*10^-6}]]

sol00 = NDSolve[{X''[NN] - (FN[NN, 0, 10^-8]) == 0, X[0] == 0, X'[0] == 0}, X, {NN, 0, 140}]

Plot[X[NN] /. {sol00}, {NN, 0, 10}, Evaluated -> True]
Improve this question
$\endgroup$
4
  • $\begingroup$ It could be better you show your function FN without this code. $\endgroup$
    – Alex Trounev
    Oct 17, 2020 at 16:18
  • $\begingroup$ @AlexTrounev Yes, I tried to do so, but the function is pretty big and ugly! It would have been a mess and perhaps less clear. $\endgroup$
    – user1886681
    Oct 20, 2020 at 18:24
  • $\begingroup$ What do you get after day or two of waiting? $\endgroup$
    – Alex Trounev
    Oct 21, 2020 at 16:20
  • $\begingroup$ @AlexTrounev I get similar results, very nice approximation! $\endgroup$
    – user1886681
    Oct 23, 2020 at 21:51

1 Answer 1

Reset to default
1
+50
$\begingroup$

We can decrease evaluation time to few minutes by filtering function FN as follows:

(*constants*)e = -1.60217733*10^-19;
m = 9.109389699999999*10^-31;
epsilon = 8.854187817620391*10^-12; lw = 0.026; kk = 1; beta = 
 Sqrt[1 - 1/(4000/0.511)^2]; gamma = 4000/0.511; c = 
 3*10^8; kw = (2 Pi)/0.026; sigma = 
 10^(-8); coeff = ((e*lw^2)/(gamma*m*beta^2*c^2))*
      (1/(10^10*e)/((2*Pi*(30/10^6))/10^5))*
      Exp[-(sigma^2/(2*(10^(-5))^2))]; 

(*basic equations*)

rs2 = {zp, xp + kk/(gamma*kw) Sin[kw*zp], 0};
ro2 = {(nn + 10000)*lw, x + kk/(gamma*kw) Sin[kw*(nn + 10000)*lw], 0};

betas = {beta - kk^2/(4 gamma^2) Cos[2 kw*zp], kk/gamma Cos[kw*zp], 0};
betao = {beta - kk^2/(4 gamma^2) Cos[2 kw*(nn + 10000)*lw], 
   kk/gamma Cos[kw*(nn + 10000)*lw], 0};

betaDot = {(c*kk^2*kw)/(2 gamma^2) Sin[
     2 kw*zp], -((c*kk*kw)/gamma) Sin[kw*zp], 0};

deltar2 = ro2 - rs2;
Rgam2 = Sqrt[deltar2[[1]]^2 + deltar2[[2]]^2];
Ec2 = (e/(4 Pi*epsilon)) (deltar2/Rgam2 - 
      betas)/(gamma^2 Rgam2^2 (1 - (deltar2/Rgam2).betas)^3);
Erad2 = (e/(4 Pi*epsilon)) Cross[deltar2/Rgam2, 
     Cross[deltar2/Rgam2 - betas, betaDot]]/(c*
      Rgam2*(1 - (deltar2/Rgam2).betas)^3);

Bc2 = Cross[deltar2/Rgam2, Ec2];
Brad2 = Cross[deltar2/Rgam2, Erad2];
Fbc2 = Cross[betao, Bc2];
Fbrad2 = Cross[betao, Brad2];
sumEtran = (Ec2[[2]] + Erad2[[2]]);
sumFBtran = Fbc2[[2]] + Fbrad2[[2]];
ZPRIME[nn_?NumericQ, x_?NumericQ] := 
    zp /. FindRoot[sigma == (1/(gamma*kw))*Sqrt[gamma^2 + kk^2]*
              (EllipticE[kw*(nn + 10000)*lw, kk^2/(gamma^2 + kk^2)] - 
                 EllipticE[kw*zp, kk^2/(gamma^2 + kk^2)]) - 
            beta*Sqrt[((nn + 10000)*lw - zp)^2 + 
                  (x + (kk*Sin[kw*(nn + 10000)*lw])/(gamma*kw) - 
                       (kk*Sin[kw*zp])/(gamma*kw))^2], {zp, 0}];

FNz = coeff*(sumEtran + sumFBtran) /. 
     {zp -> ZPRIME[nn, x-xp]};

Now instead of

FN[n_?NumericQ] :=   
     NIntegrate[
   Exp[-(xp^2/(2*(30/10^6)^2))]*
    Evaluate[FNz /. {x -> 0, xp -> xp, nn -> n}], 
       {xp, -300/10^6, 300/10^6}];

we use filtered function fp based on list interpolation. First we recognize that function fp is periodic with a period of 1

lst1 = Table[{n, 
   NIntegrate[
     Exp[-(xp^2/(2*(30/10^6)^2))]*
      Evaluate[FNz /. {x -> 0, xp -> xp, nn -> n}], 
         {xp, -300/10^6, 300/10^6}, PrecisionGoal -> 5] // Quiet}, {n, 0, 1,.005}];
lst2 = Table[{n, 
       NIntegrate[
         Exp[-(xp^2/(2*(30/10^6)^2))]*
          Evaluate[FNz /. {x -> 0, xp -> xp, nn -> n}], 
             {xp, -300/10^6, 300/10^6}, PrecisionGoal -> 5] // Quiet}, {n, 1,3,.02}];
        ListPlot[{lst1,lst2}]

Figure 1 So we can make periodic interpolation as follows

fp = Interpolation[Join[lst1, {{1, lst1[[1, 2]]}}], 
  PeriodicInterpolation -> True]

With this function we integrate equation as

sol00 = NDSolve[{X''[n] - fp[n] == 0, 
       X[0] == 0, X'[0] == 0}, X, {n, 0, 140}]

Plot[X[nn] /. {sol00}, {nn, 0, 140},Frame -> True, FrameLabel -> {"N", "X"}]

Figure 2

Finally we can test how periodic interpolation is good for this problem. We calculate 160 points in the beginning and 60 random points in the end of interval {NN,0,160}, and compare points with fp. We can check that only 3 points from 220 not follow to fp. Therefore it is good approximation.
Figure 3

Improve this answer
$\endgroup$
9
  • $\begingroup$ Jugging from OPs question I reckon that ZPRIME depends on xprime (the integration variable in F) maybe OP can comment on this issue with this running example at hand. Further I would strongly recommend circumventing the intermediate symbolic expressions by introducing numeric functions. One should also check that the computed integrals are accurate (enough) a PrecisionGoal -> 5 and Quiet might kill features and possibly critical error messages. What this solution boils down to with the current parameters is X(N)=(c N^2)/2, where c is the average of F(N) over the computational domain. $\endgroup$
    – N0va
    Oct 20, 2020 at 23:14
  • $\begingroup$ Rerunning the code with more gird points (e.g. lst =Table[{n,...}, {n, 0, 140,1/2}];) for the interpolation completely changes the result, since F is highly oscillatory. OP should check the given expressions and parameters and I would recommend getting F under control before attempting a solution of the ODE.If the displayed behavior of F is correct I would recommend trying to find a good value for its average (c) over the computational domain and approximate it as constant and use the analytical solution X(N)=(c N^2)/2 in this approximation. $\endgroup$
    – N0va
    Oct 20, 2020 at 23:34
  • $\begingroup$ @N0va Yes, ZPRIME is a function of xprime. I am checking the code/example right now. The xprime or "xp" as you have used Alex is missing from the ZPRIME[nn_?NumericQ, x_?NumericQ] function. Check my ZPRIME function in my original post. $\endgroup$
    – user1886681
    Oct 21, 2020 at 0:06
  • $\begingroup$ @user1886681 No, xp is not missing, since ZPRIME depends on x-xp. We then can put FNz = coeff*(sumEtran + sumFBtran) /. {zp -> ZPRIME[nn, x-xp]};. And so on. $\endgroup$
    – Alex Trounev
    Oct 21, 2020 at 0:14
  • $\begingroup$ @N0va I have calculated FN with zp -> ZPRIME[nn, x-xp] and result for X. You are right that it is not stable. Probably we should accumulate FN with different set of points. $\endgroup$
    – Alex Trounev
    Oct 21, 2020 at 10:24

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.