Solution numérique d'une équation intégro-différentielle

Sep 18 2020

Je voudrais résoudre numériquement l'équation intégro-différentielle suivante $$ \partial_t \rho(t,x) \,=\, \partial_x\big(f'(x)\,\rho(t,x)\big) \int_0^\infty f(\xi)\,\rho(t,\xi)\,d\xi \;+\\ +\; \partial_x\big(g'(x)\,\rho(t,x)\big) \int_0^\infty g(\xi)\,\rho(t,\xi)\,d\xi $$ où:

$\rho$ est une distribution de probabilité sur $[0,\infty)$ qui peut en fait dégénérer en une combinaison convexe d'un delta de Dirac et d'une fonction de densité;
la condition initiale $\rho(0,x)$ peut être choisi de manière appropriée, de sorte que $\int_0^\infty\rho(0,x)\,dx=1$;
disons les fonctions $f,g$sont donnés. Ils sont strictement croissants, lisses mais pas analytiques à$0$ En effet $f^{(k)}(0)=g^{(k)}(0)=0$ pour tous $k\geq1$.

J'ai essayé avec DSolve, mais aucune solution exacte n'est trouvée. Ensuite, j'ai essayé avec NDSolve et j'obtiens l'erreur suivante:

NDSolve :: delpde: les équations différentielles partielles de retard ne sont actuellement pas prises en charge par NDSolve.

Est-il possible de résoudre cette équation en utilisant Mathematica? J'utilise Mathematica 11.

Éditer

Voici la définition de $f,g$. Laisser$L(x)$ être une fonction linéaire par morceaux prenant valeur $l_0$ pour $x\leq x_0$, $l_0+\frac{x-x_0}{x_1-x_0}\,(x_1-x_0)$ pour $x_0\leq x\leq x_1$ et $l_1$ pour $x\geq x_1$. Puis réglez:$$ E(x) = \int_{-\infty}^{\infty} L(xz)\, \frac{e^{-\frac{z^2}{2}}}{\sqrt{2\pi}}\, dz $$ enfin réparer $c$ positif, $\epsilon\in(0,1)$ et laissez $$ f(x) = c\,E\big((1+\epsilon)\,x\big)-c \quad,\quad g(x) = c\,E\big((1-\epsilon)\,x\big)+c\;. $$ Par exemple, réparer $l_0=-2.5,\,l_1=7.5,\,x_0=0.5,\,x_1=1.5$ et $c=1,\,\epsilon=0.6\,$.

Modifier 2

J'ai obtenu un tracé de la solution mettant en œuvre la méthode numérique des lignes suggérée par @bbgodfrey, mais il y a les mêmes problèmes pour $x$ proche de $0$. Voici le tracé résultant, de deux points de vue:

Solution $\rho(t,r)$obtenu par la méthode numérique des lignes. Vue 1

Solution $\rho(t,r)$obtenu par la méthode numérique des lignes. Vue 2

Il semble que quelque chose se passe autour $t\approx0.5$. Quelles sont ces lignes droites? Existe-t-il un moyen de voir clairement l'apparence d'une fonction Delta et de la distinguer des problèmes numériques?

Voici mon code:

n = 1000; rmax = 5; T = 2;
X = Table[rmax/n*(i - 1), {i, 1, n + 1}];
Rho[t_] := Table[Subscript[ρ, i][t], {i, 1, n + 1}];
F = Table[f[X[[i]] $MachinePrecision], {i, 1, n + 1}]; G = Table[g[X[[i]] $MachinePrecision], {i, 1, n + 1}];
DF = Table[Df[X[[i]] $MachinePrecision], {i, 1, n + 1}]; DG = Table[Dg[X[[i]] $MachinePrecision], {i, 1, n + 1}];

(* Initial condition *)
gamma[r_] := 1/(Gamma[k] θ^k) r^(k - 1) Exp[-r/θ]
k = 10; θ = 0.1;
ic = Thread[ Drop[Rho[0], -1] == Table[gamma[X[[i]]], {i, 1, n}] ];

(* Boundary condition *)
Subscript[ρ, n + 1][t_] := 0

(* ODE's *)
rhs[t_] := 
  ListCorrelate[{-1, 1}, DF*Rho[t]]*Total[F*Rho[t]] + 
  ListCorrelate[{-1, 1}, DG*Rho[t]]*Total[G*Rho[t]]
lhs[t_] := Drop[D[Rho[t], t] , -1]
eqns = Thread[lhs[t] == rhs[t]];

lines = 
  NDSolve[
    {eqns, ic}, Drop[Rho[t], -1], {t, 0, T}, 
    Method -> {"EquationSimplification" -> "Residual"}];

ParametricPlot3D[
  Evaluate[Table[{rmax/n*i, t, First[Subscript[ρ, i][t] /. lines]}, {i, 1, n/2}]], 
  {t, 0, 1}, 
  AxesLabel -> {"r", "t", "ρ"}, BoxRatios -> {1, 1, 1}]

Réponses

3 AlexTrounev Sep 28 2020 at 23:26

Puisque dans le code original il y a des instabilités dues à une approximation d'ordre faible, nous pouvons utiliser l'algorithme numérique de 4ème ordre que j'ai développé pour le modèle démographique Lotka-McKendrick (voir le tout dernier code dans ma réponse). Nous définissons d'abord la fonction en f, gutilisant l'expression exacte suivante pour$E(x)$:

l0 = -25/10; l1 = 75/10; x0 = 1/2; x1 = 3/2; c = 1; eps = 3/5; 
L[x_] := Piecewise[{{l0, x <= x0}, {l0 + (l1 - l0) (x - x0)/(x1 - x0),
     x0 < x <= x1}, {l1, x > x1}}]; 
Integrate[L[x z] Exp[-z^2/2], {z, -Infinity, Infinity}, 
  Assumptions -> {x > 0}]/Sqrt[2 Pi]

(*1/(4 Sqrt[2 \[Pi]])5 \[ExponentialE]^(-(9/(8 x^2))) (-\
\[ExponentialE]^((9/(8 x^2))) Sqrt[2 \[Pi]]-8 x+8 \
\[ExponentialE]^(1/x^2) x+2 \[ExponentialE]^(9/(8 x^2)) Sqrt[2 \[Pi]] \
Erf[1/(2 Sqrt[2] x)]-3 \[ExponentialE]^(9/(8 x^2)) Sqrt[2 \[Pi]] \
Erf[3/(2 Sqrt[2] x)]+3 \[ExponentialE]^(9/(8 x^2)) Sqrt[2 \[Pi]] \
Erfc[3/(2 Sqrt[2] x)])*)

Par conséquent, nous pouvons définir explicitement les fonctions $f(x),g(x),E(x),E'(x)f'(x), g'(x)$comme f,g,eL,eL1,df,dg, nous avons

eL[x_] := 
 1/(4 Sqrt[2 \[Pi]])
   5 E^(-(9/(
   8 x^2))) (-E^((9/(8 x^2))) Sqrt[2 \[Pi]] - 8 x + 8 E^(1/x^2) x + 
    2 E^(9/(8 x^2)) Sqrt[2 \[Pi]] Erf[1/(2 Sqrt[2] x)] - 
    3 E^(9/(8 x^2)) Sqrt[2 \[Pi]] Erf[3/(2 Sqrt[2] x)] + 
    3 E^(9/(8 x^2)) Sqrt[2 \[Pi]] Erfc[3/(2 Sqrt[2] x)]); 
eL1[x_] := (
  45 E^(-(9/(
    8 x^2))) (-E^((9/(8 x^2))) Sqrt[2 \[Pi]] - 8 x + 8 E^(1/x^2) x + 
     2 E^(9/(8 x^2)) Sqrt[2 \[Pi]] Erf[1/(2 Sqrt[2] x)] - 
     3 E^(9/(8 x^2)) Sqrt[2 \[Pi]] Erf[3/(2 Sqrt[2] x)] + 
     3 E^(9/(8 x^2)) Sqrt[2 \[Pi]] Erfc[3/(2 Sqrt[2] x)]))/(
  16 Sqrt[2 \[Pi]] x^3) + (
  5 E^(-(9/(
    8 x^2))) (-8 + 8 E^(1/x^2) + (9 E^(9/(8 x^2)) Sqrt[\[Pi]/2])/(
     2 x^3) + 18/x^2 - (18 E^(1/x^2))/x^2 - (
     9 E^(9/(8 x^2)) Sqrt[\[Pi]/2] Erf[1/(2 Sqrt[2] x)])/x^3 + (
     27 E^(9/(8 x^2)) Sqrt[\[Pi]/2] Erf[3/(2 Sqrt[2] x)])/(2 x^3) - (
     27 E^(9/(8 x^2)) Sqrt[\[Pi]/2] Erfc[3/(2 Sqrt[2] x)])/(2 x^3)))/(
  4 Sqrt[2 \[Pi]]); f[x_] := c eL[(1 + eps) x] - c; 
df[x_] := c (1 + eps) eL1[(1 + eps) x]; 
g[x_] := c eL[(1 - eps) x] + c; 
dg[x_] := c (1 - eps) eL1[(1 - eps) x];

Deuxième étape, nous appelons

Needs["DifferentialEquations`NDSolveProblems`"];
Needs["DifferentialEquations`NDSolveUtilities`"]; 
Get["NumericalDifferentialEquationAnalysis`"];

Maintenant, nous définissons la grille et les poids pour l'intégration numérique en utilisant GaussianQuadratureWeights[]et DifferentiationMatrixsur la même grille en utilisant FiniteDifferenceDerivative:

np = 100; gqw = GaussianQuadratureWeights[np, 0, 5];
ugrid = gqw[[All, 1]]; weights = gqw[[All, 2]]; fd = 
 NDSolve`FiniteDifferenceDerivative[Derivative[1], ugrid]; m = 
 fd["DifferentiationMatrix"];

Enfin, nous définissons tous les vecteurs, matrices, équations nécessaires et résolvons le système d'ODE en utilisant NDSolve

Quiet[varf = Table[df[ugrid[[i]]] u[i][t], {i, Length[ugrid]}]; 
 varg = Table[dg[ugrid[[i]]] u[i][t], {i, Length[ugrid]}]; 
 varu = Table[u[i][t], {i, Length[ugrid]}]; 
 var = Table[u[i], {i, Length[ugrid]}]; ufx = m.varf; ugx = m.varg; 
 intf = Table[f[ugrid[[i]]] weights[[i]], {i, np}]; 
 intg = Table[g[ugrid[[i]]] weights[[i]], {i, np}]]; 
u0[r_] := 1/(Gamma[k] \[Theta]^k) r^(k - 1) Exp[-r/\[Theta]]
k = 10; \[Theta] = 0.1;

ics = Table[u[i][0] == u0[ugrid[[i]]], {i, np}]; eqns = 
 Table[D[u[i][t], t] == 
   ufx[[i]] (intf.varu) + ugx[[i]] (intg.varu), {i, np}]; tmax = 2;
sol = NDSolve[{eqns, ics}, var, {t, 0, tmax}, 
   Method -> {"EquationSimplification" -> "Residual"}];

Visualisation de la solution numérique

lst = Flatten[
   Table[{t, ugrid[[i]], u[i][t] /. sol[[1]]}, {t, 0, 2, 1/50}, {i, 
     np}], 1];
ListPlot3D[lst, Mesh -> None, PlotRange -> All, 
 AxesLabel -> {"t", "x"}]

Nous pouvons comparer ce résultat avec le code original fonctionnant pour n=50(image de gauche) et n=100(droite). Sur l'image de gauche, nous pouvons reconnaître la solution indiquée ci-dessus. Mais il y a aussi des oscillations non physiques avec une amplitude augmentant 10 fois avec n augmente de 50 à 100. Code original tel que je suis utilisé pourn=50

eL[x_] := 
 1/(4 Sqrt[2 \[Pi]])
   5 E^(-(9/(
   8 x^2))) (-E^((9/(8 x^2))) Sqrt[2 \[Pi]] - 8 x + 8 E^(1/x^2) x + 
    2 E^(9/(8 x^2)) Sqrt[2 \[Pi]] Erf[1/(2 Sqrt[2] x)] - 
    3 E^(9/(8 x^2)) Sqrt[2 \[Pi]] Erf[3/(2 Sqrt[2] x)] + 
    3 E^(9/(8 x^2)) Sqrt[2 \[Pi]] Erfc[3/(2 Sqrt[2] x)]); 
eL1[x_] := (
  45 E^(-(9/(
    8 x^2))) (-E^((9/(8 x^2))) Sqrt[2 \[Pi]] - 8 x + 8 E^(1/x^2) x + 
     2 E^(9/(8 x^2)) Sqrt[2 \[Pi]] Erf[1/(2 Sqrt[2] x)] - 
     3 E^(9/(8 x^2)) Sqrt[2 \[Pi]] Erf[3/(2 Sqrt[2] x)] + 
     3 E^(9/(8 x^2)) Sqrt[2 \[Pi]] Erfc[3/(2 Sqrt[2] x)]))/(
  16 Sqrt[2 \[Pi]] x^3) + (
  5 E^(-(9/(
    8 x^2))) (-8 + 8 E^(1/x^2) + (9 E^(9/(8 x^2)) Sqrt[\[Pi]/2])/(
     2 x^3) + 18/x^2 - (18 E^(1/x^2))/x^2 - (
     9 E^(9/(8 x^2)) Sqrt[\[Pi]/2] Erf[1/(2 Sqrt[2] x)])/x^3 + (
     27 E^(9/(8 x^2)) Sqrt[\[Pi]/2] Erf[3/(2 Sqrt[2] x)])/(2 x^3) - (
     27 E^(9/(8 x^2)) Sqrt[\[Pi]/2] Erfc[3/(2 Sqrt[2] x)])/(2 x^3)))/(
  4 Sqrt[2 \[Pi]]); f[x_] := c eL[(1 + eps) x] - c; 
df[x_] := c (1 + eps) eL1[(1 + eps) x]; 
g[x_] := c eL[(1 - eps) x] + c; dg[x_] := c (1 - eps) eL1[(1 - eps) x];

n = 50; rmax = 5; T = 2;
X = Table[rmax/n*(i - 1) + 10^-6, {i, 1, n + 1}];
Rho[t_] := Table[Subscript[\[Rho], i][t], {i, 1, n + 1}];
F = Table[f[X[[i]] ], {i, 1, n + 1}];
G = Table[g[X[[i]] ], {i, 1, n + 1}];
DF = Table[df[X[[i]]], {i, 1, n + 1}];
DG = Table[dg[X[[i]] ], {i, 1, n + 1}];

(*Initial condition*)
gamma[r_] := 1/(Gamma[k] \[Theta]^k) r^(k - 1) Exp[-r/\[Theta]]
k = 10; \[Theta] = 0.1;
ic = Thread[Drop[Rho[0], -1] == Table[gamma[X[[i]]], {i, 1, n}]];

(*Boundary condition*)
Subscript[\[Rho], n + 1][t_] := 0

(*ODE's*)
rhs[t_] := 
 ListCorrelate[{-1, 1}, DF*Rho[t]]*Total[F*Rho[t]] + 
  ListCorrelate[{-1, 1}, DG*Rho[t]]*Total[G*Rho[t]]
lhs[t_] := Drop[D[Rho[t], t], -1]
eqns = Thread[lhs[t] == rhs[t]];

lines = NDSolve[{eqns, ic}, Drop[Rho[t], -1], {t, 0, T}, 
   Method -> {"EquationSimplification" -> "Residual"}];

Visualisation des solutions numériques pour n=50(gauche) et n=100(droite)

lst = Table[{t, X[[i]], Subscript[\[Rho], i][t] /. lines[[1]]}, {t, 0,
     T, 1/25}, {i, n}];

ListPlot3D[Flatten[lst, 1], ColorFunction -> "Rainbow", Mesh -> None, 
 AxesLabel -> {"t", "x", ""}, PlotRange -> All]

1 Dominic Sep 23 2020 at 20:16

Ce n'est pas une réponse, mais quelques commentaires sur la résolution de ce type de problème qui sont trop longs et à faire en commentaire à la question.

Concernant la mise à l'échelle vers le haut et vers le bas: à mon avis, pour devenir compétent dans la résolution de problèmes difficiles, il est impératif d'apprendre à réduire le problème, puis à le remonter. Par exemple, vous avez:$$ \frac{\partial \rho}{\partial t}=\frac{\partial}{\partial t}\left(f'\rho\right)\int_0^{\infty} f(x)\rho(t,x)dx+\cdots $$Remarquez les points. Une fois supprimé, cela le réduit à une forme plus simple. Pouvez-vous résoudre celui-là? Peut-être qu'il n'a pas de solution. Je ne sais pas. Que diriez-vous de retirer le$f'\rho$ terme, dites:

$$ \frac{\partial \rho}{\partial t}+\frac{\partial p}{\partial x}=\int_0^{\infty} f(x)\rho(t,x)dx $$

Celui-là? Que diriez-vous de retirer le$f(x)$terme dans l'intégrale? Combien devez-vous réduire tout en conservant sa nature PIDE afin de le résoudre? Que diriez-vous de résoudre n'importe quel PIDE simple (quelque peu similaire) pour perfectionner la méthode, puis d'ajouter de la complexité (des termes) au problème jusqu'à ce que vous atteigniez l'équation que vous souhaitez résoudre.

Bien sûr, cela demande beaucoup de travail et parfois, vous rencontrerez bien sûr des problèmes où son augmentation supplémentaire entraîne une difficulté importante à résoudre. Mais étonnamment, cette méthode a souvent été très efficace avec des problèmes difficiles sur lesquels j'ai travaillé, mais pas toujours. Voici un exemple:$$ f+\frac{\partial f}{\partial x}+\frac{\partial f}{\partial y}=\int_x^{\infty} \int_y^{\infty}f(u,v)dudv $$beau hein, mais un peu intimidant. Que diriez-vous de le réduire:$$ f+\frac{df}{dx}=\int_x^{\infty} f(u)du $$C'est plus facile et il s'avère que la solution à celle-là mène facilement à la solution de la première. :)