# NDSolve with embedded FindRoot and NIntegrate

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[]^2 + deltar2[]^2];
Ec2 = (e/(4 Pi*epsilon)) (deltar2/Rgam2 - betas)/(gamma^2 Rgam2^2 (1 - (deltar2/Rgam2).betas)^3);

Bc2 = Cross[deltar2/Rgam2, Ec2];

Fbc2 = Cross[betao, Bc2];

(*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, X' == 0}, X, {NN, 0, 140}]

Plot[X[NN] /. {sol00}, {NN, 0, 10}, Evaluated -> True]

• It could be better you show your function FN without this code. Oct 17, 2020 at 16:18
• @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. Oct 20, 2020 at 18:24
• What do you get after day or two of waiting? Oct 21, 2020 at 16:20
• @AlexTrounev I get similar results, very nice approximation! Oct 23, 2020 at 21:51

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};

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

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

Bc2 = Cross[deltar2/Rgam2, Ec2];
Fbc2 = Cross[betao, Bc2];
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]};


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}] 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, X' == 0}, X, {n, 0, 140}]

Plot[X[nn] /. {sol00}, {nn, 0, 140},Frame -> True, FrameLabel -> {"N", "X"}] 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. • 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.
– N0va
Oct 20, 2020 at 23:14
• 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.
– N0va
Oct 20, 2020 at 23:34
• @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. Oct 21, 2020 at 0:06
• @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. Oct 21, 2020 at 0:14
• @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. Oct 21, 2020 at 10:24