NDSolve avec FindRoot et NIntegrate intégrés
J'essaie de tracer des résultats pour une équation différentielle très simple de la forme:
$$\frac{\partial^2 x(N,z'(N))}{\partial N^2} = F(N,z'(N)), $$
où $z'(N)$ est une fonction de $N$ qui doit être résolu en utilisant FindRoot pour à chaque $N$ position, et $F(N,z')$ est une mauvaise équation qui résulte d'une intégration numérique sur:
$$ F(N,z') = \int_{-\infty}^{\infty} \exp\left( -\frac{x'^2}{2\sigma_{x'}^2} \right) F(N,z',x')dx'$$.
Donc, j'ai assemblé du code mathématique mais il tourne très lentement (de l'ordre d'un jour ou deux)! J'ai remarqué que certaines choses affectaient la vitesse du code, en particulier le coefficient numérique devant$F(N,z'(N))$. Mais je me demandais s'il y avait de l'aide à apporter pour obtenir des résultats meilleurs / plus rapides! Toute aide serait grandement appréciée!
Remarque: j'ai dû utiliser $NN$ au lieu de $N$dans mes équations parce que c'est une fonction en mathématique. De plus, dans la fonction FN, je dois copier et coller la sortie de FNzprime (un vilain gâchis) dans l'intégrande pour qu'elle soit évaluée.
(*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]
Réponses
Nous pouvons réduire le temps d'évaluation à quelques minutes en filtrant la fonction FNcomme suit:
(*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]};
Maintenant au lieu de
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}];
nous utilisons une fonction filtrée fpbasée sur l'interpolation de liste. Tout d'abord, nous reconnaissons que la fonction fpest périodique avec une période de 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}]
fp = Interpolation[Join[lst1, {{1, lst1[[1, 2]]}}],
PeriodicInterpolation -> True]
Avec cette fonction, nous intégrons l'équation comme
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"}]
Enfin, nous pouvons tester comment l'interpolation périodique est bonne pour ce problème. Nous calculons 160 points au début et 60 points aléatoires à la fin de l'intervalle {NN,0,160}, et comparons les points avec fp. Nous pouvons vérifier que seuls 3 points de 220 ne suivent pas fp. C'est donc une bonne approximation.
