NDSolve mit eingebettetem FindRoot und NIntegrate

Question

Ich versuche, Ergebnisse für eine sehr einfache Differentialgleichung der Form zu zeichnen:

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

wo $z'(N)$ ist eine Funktion von $N$ das muss mit FindRoot für jeden gelöst werden $N$ Position und $F(N,z')$ ist eine böse Gleichung, die sich aus einer numerischen Integration ergibt über:

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

Also habe ich etwas Mathematica-Code zusammengestellt, aber er läuft furchtbar langsam (in der Größenordnung von ein oder zwei Tagen)! Ich bemerkte, dass es einige Dinge gab, die die Geschwindigkeit des Codes beeinflussten, insbesondere den numerischen Koeffizienten vor$F(N,z'(N))$. Aber ich habe mich gefragt, ob es Hilfe gibt, um bessere / schnellere Ergebnisse zu erzielen! Jede Hilfe wäre sehr dankbar!

Hinweis: Ich musste verwenden $NN$ anstelle von $N$in meinen Gleichungen, weil es eine Funktion in der Mathematik ist. Außerdem muss ich in der FN-Funktion die Ausgabe von FNzprime (ein hässliches Durcheinander) kopieren und in den Integranden einfügen, damit sie ausgewertet wird.

(*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]

Answer 1

Wir können die Auswertungszeit auf einige Minuten verkürzen, indem wir die Funktion FNwie folgt filtern :

(*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]};

Jetzt statt

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

Wir verwenden eine gefilterte Funktion fpbasierend auf Listeninterpolation. Zuerst erkennen wir, dass die Funktion fpmit einer Periode von 1 periodisch ist

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

Wir können also die periodische Interpolation wie folgt durchführen

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

Mit dieser Funktion integrieren wir die Gleichung als

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"}]

Schließlich können wir testen, wie gut die periodische Interpolation für dieses Problem ist. Wir berechnen 160 Punkte am Anfang und 60 zufällige Punkte am Ende des Intervalls {NN,0,160}und vergleichen die Punkte mit fp. Wir können überprüfen, dass nur 3 Punkte von 220 nicht folgen fp. Daher ist es eine gute Annäherung.