Trois PDE couplés à résoudre de manière semi-analytique / analytique
J'ai essayé de résoudre les trois PDE couplés suivants où l'objectif final est de trouver les distributions $\theta_h, \theta_c$ et $\theta_w$:
$x\in[0,1]$ et $y\in[0,1]$
$$\frac{\partial \theta_h}{\partial x}+\beta_h (\theta_h-\theta_w) = 0 \tag A$$
$$\frac{\partial \theta_c}{\partial y} + \beta_c (\theta_c-\theta_w) = 0 \tag B$$
$$\lambda_h \frac{\partial^2 \theta_w}{\partial x^2} + \lambda_c V\frac{\partial^2 \theta_w}{\partial y^2}-\frac{\partial \theta_h}{\partial x} - V\frac{\partial \theta_c}{\partial y} = 0 \tag C$$
où, $\beta_h, \beta_c, V, \lambda_h, \lambda_c$sont des constantes. Les conditions aux limites sont:
$$\frac{\partial \theta_w(0,y)}{\partial x}=\frac{\partial \theta_w(1,y)}{\partial x}=\frac{\partial \theta_w(x,0)}{\partial y}=\frac{\partial \theta_w(x,1)}{\partial y}=0$$
$$\theta_h(0,y)=1, \theta_c(x,0)=0$$
Un utilisateur d'échange de pile de mathématiques m'a suggéré les étapes suivantes qui pourraient contribuer à résoudre ce problème:
- Représenter chacune des trois fonctions à l'aide de séries de Fourier 2D
- Observez que toutes les équations sont linéaires
- Il n'y a donc pas de couplage de fréquence
- Ainsi pour chaque paire de fréquences $\omega_x$, $\omega_y$ il y aura une solution à partir d'une combinaison linéaire de seulement ces termes
- Appliquer les conditions aux limites directement à chacune des trois séries
- Notez que par orthogonalité la condition aux limites doit s'appliquer à chaque terme de la série de Fourier
- Branchez la série de Fourier dans PDE et résolvez la correspondance des coefficients ( voir ici par exemple en 1D ). Assurez-vous de traiter séparément les cas où l'une ou les deux fréquences sont nulles.
- Si vous considérez toutes les équations pour une paire de fréquences donnée, vous pouvez les organiser en une équation $M\alpha = 0$, où $\alpha$ sont des coefficients de Fourier pour ces fréquences, et $M$ est une petite matrice clairsemée (qch comme 12x12) qui ne dépendra que des constantes.
- Pour chaque fréquence, les solutions autorisées seront dans l'espace nul de cette matrice. Si vous ne parvenez pas à résoudre analytiquement l'espace nul, ce n'est pas un gros problème - le calcul numérique de l'espace nul est facile, en particulier pour les petites matrices.
Quelqu'un peut-il m'aider à appliquer ces étapes dans Mathematica?
PDE1 = D[θh[x, y], x] + bh*(θh[x, y] - θw[x, y]) == 0;
PDE2 = D[θc[x, y], y] + bc*(θc[x, y] - θw[x, y]) == 0;
PDE3 = λh*D[θw[x, y], {x, 2}] + λc*V*(D[θw[x, y], {y, 2}]) - D[θh[x, y], x] - V*D[θc[x, y], y] ==0
bh=0.433;bc=0.433;λh = 2.33 10^-6; λc = 2.33 10^-6; V = 1;
Solution NDSolve (mauvais résultats)
PDE1 = D[θh[x, y], x] + bh*(θh[x, y] - θw[x, y]) == 0;
PDE2 = D[θc[x, y], y] + bc*(θc[x, y] - θw[x, y]) == 0;
PDE3 = λh*D[θw[x, y], {x, 2}] + λc*V*(D[θw[x, y], {y, 2}]) - D[θh[x, y], x] - V*D[θc[x, y], y] == NeumannValue[0, x == 0.] + NeumannValue[0, x == 1] +
NeumannValue[0, y == 0] + NeumannValue[0, y == 1];
bh = 1; bc = 1; λh = 1; λc = 1; V = 1;(*Random \
values*)
sol = NDSolve[{PDE1, PDE2, PDE3, DirichletCondition[θh[x, y] == 1, x == 0], DirichletCondition[θc[x, y] == 0, y == 0]}, {θh, θc, θw}, {x, 0, 1}, {y, 0, 1}]
Plot3D[θw[x, y] /. sol, {x, 0, 1}, {y, 0, 1}]
Vers une solution séparable
J'ai écrit $\theta_h(x,y) = \beta_h e^{-\beta_h x} \int e^{\beta_h x} \theta_w(x,y) \, \mathrm{d}x$ et $\theta_c(x,y) = \beta_c e^{-\beta_c y} \int e^{\beta_c y} \theta_w(x,y) \, \mathrm{d}y$ et éliminé $\theta_h$ et $\theta_c$de Eq. (C). Puis j'ai utilisé l'ansatz$\theta_w(x,y) = e^{-\beta_h x} f(x) e^{-\beta_c y} g(y)$sur ce nouvel Eq. (C) pour le séparer en$x$ et $y$Composants. Puis en utilisant$F(x) := \int f(x) \, \mathrm{d}x$ et $G(y) := \int g(y) \, \mathrm{d}y$, J'obtiens les deux équations suivantes:
\ begin {eqnarray} \ lambda_h F '' '- 2 \ lambda_h \ beta_h F' '+ \ left ((\ lambda_h \ beta_h - 1) \ beta_h - \ mu \ right) F' + \ beta_h ^ 2 F & = & 0, \\ V \ lambda_c G '' '- 2 V \ lambda_c \ beta_c G' '+ \ gauche ((\ lambda_c \ beta_c - 1) V \ beta_c + \ mu \ right) G' + V \ beta_c ^ 2 G & = & 0, \ end {eqnarray} avec une constante de séparation$\mu \in \mathbb{R}$. Cependant, je ne pouvais pas aller plus loin.
Une équation différentielle intégrale partielle
Éliminer $\theta_h, \theta_c$de l'Eq. (C) donne lieu à une équation différentielle partio-intégrale:
\ begin {eqnarray} 0 & = & e ^ {- \ beta_h x} \ left (\ lambda_h e ^ {\ beta_h x} \ frac {\ partial ^ 2 \ theta_w} {\ partial x ^ 2} - \ beta_h e ^ {\ beta_h x} \ theta_w + \ beta_h ^ 2 \ int e ^ {\ beta_h x} \ theta_w \, \ mathrm {d} x \ right) + \\ && + V e ^ {- \ beta_c y} \ gauche (\ lambda_c e ^ {\ beta_c y} \ frac {\ partial ^ 2 \ theta_w} {\ partial y ^ 2} - \ beta_c e ^ {\ beta_c y} \ theta_w + \ beta_c ^ 2 \ int e ^ { \ beta_c y} \ theta_w \, \ mathrm {d} y \ right). \ end {eqnarray}
POINTS
Pour bc = 4; bh = 2; λc = 0.01; λh = 0.01; V = 2;
Cependant, les mêmes paramètres mais avec V=1
fonctionnent bien.

Du matériel de référence pour les futurs utilisateurs
Afin de comprendre l'évaluation des coefficients de Fourier en utilisant le concept de minimisation des moindres carrés que @bbgodfrey utilise dans sa réponse, les futurs utilisateurs peuvent se pencher sur cet article de R. Kelman (1979). Alternativement, cette présentation et cette vidéo sont également des références utiles.
Réponses
Modifications: remplacement de l'expansion à 1 terme par une expansion à n termes; amélioration de la généralité des calculs des valeurs propres et des coefficients; code réorganisé et simplifié.
En commençant par cet ensemble d'équations, procédez comme suit pour obtenir une solution presque symbolique.
ClearAll[Evaluate[Context[] <> "*"]]
eq1 = D[θh[x, y], x] + bh (θh[x, y] - θw[x, y])
eq2 = D[θc[x, y], y] + bc (θc[x, y] - θw[x, y])
eq3 = λh D[θw[x, y], x, x] + λc V D[θw[x, y], y, y] + bh (θh[x, y] - θw[x, y]) +
V bc (θc[x, y] - θw[x, y])
Commencez par convertir ces équations en ODE par la méthode de séparation des variables.
th = Collect[(eq1 /. {θh -> Function[{x, y}, θhx[x] θhy[y]],
θw -> Function[{x, y}, θwx[x] θwy[y]]})/(θhy[y] θwx[x]),
{θhx[x], θhx'[x], θwy[y]}, Simplify];
1 == th[[1 ;; 3 ;; 2]];
eq1x = Subtract @@ Simplify[θwx[x] # & /@ %] == 0
1 == -th[[2]];
eq1y = θhy[y] # & /@ %
(* bh θhx[x] - θwx[x] + θhx'[x] == 0
θhy[y] == bh θwy[y] *)
tc = Collect[(eq2 /. {θc -> Function[{x, y}, θcx[x] θcy[y]],
θw -> Function[{x, y}, θwx[x] θwy[y]]})/(θcx[x] θwy[y]),
{θcy[y], θcy'[y], θwy[y]}, Simplify];
1 == -tc[[1]];
eq2x = θcx[x] # & /@ %
1 == tc[[2 ;; 3]];
eq2y = Subtract @@ Simplify[θwy[y] # & /@ %] == 0
(* θcx[x] == bc θwx[x]
bc θcy[y] - θwy[y] + [θcy[y] == 0 *)
tw = Plus @@ ((List @@ Expand[eq3 /. {θh -> Function[{x, y}, θhx[x] θhy[y]],
θc -> Function[{x, y}, θcx[x] θcy[y]], θw -> Function[{x, y}, θwx[x] θwy[y]]}])/
(θwx[x] θwy[y]) /. Rule @@ eq2x /. Rule @@ eq1y);
sw == -tw[[1 ;; 5 ;; 2]];
eq3x = Subtract @@ Simplify[θwx[x] # & /@ %] == 0
sw == tw[[2 ;; 6 ;; 2]];
eq3y = -Subtract @@ Simplify[θwy[y] # & /@ %] == 0
(* bh^2 θhx[x] - bh θwx[x] + sw θwx[x] + λh θwx''[x] == 0
bc^2 V θcy[y] - (sw + bc V) θwy[y] + V λc θwy''[y] == 0 *)
Avec les équations séparées en ODE, résolvez les équations dépendantes de y avec les conditions aux limites appliquées. Les expressions résultantes, impliquant RootSum
, sont longues et ne sont donc pas reproduites ici.
sy = DSolveValue[{eq2y, eq3y, θcy[0] == 0, θwy'[0] == 0}, {θwy[y], θcy[y], θwy'[1]},
{y, 0, 1}] /. C[2] -> 1;
Ceci est, bien sûr, un problème avec des solutions non triviales valeur propre que pour des valeurs discrètes de la constante de séparation, sw
. La relation de dispersion pour sw
est donnée par θwy'[1] == 0
. La x
dépendance correspondante est déterminée pour chaque valeur propre par
sx = DSolveValue[{eq1x, eq3x, θwx'[0] == 0, θwx'[1] == 0, θhx[0] == 1},
{θwx[x], θhx[x]}, {x, 0, 1}];
et c'est à ce stade que la condition aux limites non homogène,, θhx[0] == 1
est appliquée. Ce résultat est également trop long pour être reproduit ici.
Ensuite, déterminez numériquement les premières n = 6
valeurs propres (ici ), ce qui nécessite de spécifier les paramètres:
bc = 1; bh = 1; λc = 1; λh = 1; V = 1;
disp = sy[[3]]
(* RootSum[sw + #1 + sw #1 - #1^2 - #1^3 &,
(E^#1 sw + E^#1 #1 + E^#1 sw #1)/(-1 - sw + 2 #1 + 3 #1^2) &] *)
n = 6;
Plot[disp, {sw, -300, 10}, AxesLabel -> {sw, "disp"},
LabelStyle -> {15, Bold, Black}, ImageSize -> Large]

Les premières valeurs propres sont estimées à partir des zéros du tracé, puis calculées avec une grande précision.
Partition[Union @@ Cases[%, Line[z_] -> z, Infinity], 2, 1];
Reverse[Cases[%, {{z1_, z3_}, {z2_, z4_}} /; z3 z4 < 0 :> z1]][[1 ;; n]];
tsw = sw /. Table[FindRoot[disp, {sw, sw0}], {sw0, %}]
(* {-0.635232, -10.7982, -40.4541, -89.8156, -158.907, -247.736} *)
et les fonctions propres correspondantes obtenues en branchant ces valeurs de sw
into sy[1;;2]
et sx
.
Plot[Evaluate@ComplexExpand@Replace[sy[[1]],
{sw -> #} & /@ tsw, Infinity], {y, 0, 1}, AxesLabel -> {y, θwy},
LabelStyle -> {15, Bold, Black}, ImageSize -> Large]
Plot[Evaluate@ComplexExpand@Replace[sy[[2]],
{sw -> #} & /@ tsw, Infinity], {y, 0, 1}, AxesLabel -> {y, θhy},
LabelStyle -> {15, Bold, Black}, ImageSize -> Large]


Plot[Evaluate@ComplexExpand@Replace[sx[[1]],
{sw -> #} & /@ tsw, Infinity], {x, 0, 1}, AxesLabel -> {x, θwx},
LabelStyle -> {15, Bold, Black}, ImageSize -> Large, PlotRange -> {0, 1}]
Plot[Evaluate@ComplexExpand@Replace[sx[[2]],
{sw -> #} & /@ tsw, Infinity], {x, 0, 1}, AxesLabel -> {x, θhx},
LabelStyle -> {15, Bold, Black}, ImageSize -> Large, PlotRange -> {0, 1}]


Avec les premières n
fonctions propres complètes calculées, leurs coefficients sont ensuite déterminés, de sorte qu'ils peuvent être additionnés pour approcher la solution des équations d'origine. Cela se fait par les moindres carrés, car le système ODE n'est pas auto-adjoint.
syn = ComplexExpand@Replace[bh sy[[1]] /. C[2] -> 1, {sw -> #} & /@ tsw,
Infinity] // Chop//Chop;
Integrate[Expand[(1 - Array[c, n].syn)^2], {y, 0, 1}] // Chop;
coef = ArgMin[%, Array[c, n]]
(* {0.974358, 0.0243612, 0.000807808, 0.000341335, 0.0000506603, \
0,0000446734} *)
La qualité de l'ajustement est très bonne.
Plot[coef.syn - 1, {y, 0, 1}, AxesLabel -> {y, err},
LabelStyle -> {15, Bold, Black}, ImageSize -> Large]

Enfin, construisez la solution.
solw = coef.ComplexExpand@Replace[sy[[1]] sx[[1]], {sw -> #} & /@ tsw, Infinity];
Plot3D[solw, {x, 0, 1}, {y, 0, 1}, AxesLabel -> {x, y, θw},
LabelStyle -> {15, Bold, Black}, ImageSize -> Large]

solh = coef.ComplexExpand@Replace[bh sy[[1]] sx[[2]], {sw -> #} & /@ tsw, Infinity];
Plot3D[solh, {x, 0, 1}, {y, 0, 1}, AxesLabel -> {x, y, θh},
LabelStyle -> {15, Bold, Black}, ImageSize -> Large, PlotRange -> {0, 1}]

solc = coef.ComplexExpand@Replace[bc sy[[2]] sx[[1]], {sw -> #} & /@ tsw, Infinity];
Plot3D[solc, {x, 0, 1}, {y, 0, 1}, AxesLabel -> {x, y, θc},
LabelStyle -> {15, Bold, Black}, ImageSize -> Large, PlotRange -> {0, 1}]

Cette dérivation étant longue, nous montrons ici que les équations elles-mêmes sont satisfaites à l'identique.
Chop@Simplify[{eq1, eq2, eq3} /. {θh -> Function[{x, y}, Evaluate@solh],
θc -> Function[{x, y}, Evaluate@solc], θw -> Function[{x, y}, Evaluate@solw]}]
(* {0, 0, 0} *)
De plus, la condition aux limites on θh
est satisfaite à mieux que 0,004%, et la condition aux limites on θc
est satisfaite de manière identique.
Le calcul 3D correspondant a été terminé à 226346 .
La solution que j'obtiens avec la version 12.0.0 semble en effet incohérente. Je compare la solution assez proche de celle présentée sur la page de documentation pour NDSolvedans la section Problèmes possibles -> Equations différentielles partielles avec l'exemple de l'équation de Laplace avec les valeurs initiales.
Pour le système d'équation aux dérivées partielles donné et pour la valeur définie uniquement avec une seule, je peux utiliser NDSolvepour ce résultat:

La similitude n'est pas la divergence qui tombe à l'origine mais la rangée de pointes que l'on peut voir à environ $x=.3$ et $y=0.3$ pour $𝜃_h$ et $𝜃_c$. Ce couplage est bien que vraiment non physique. Mais il y a des informations apparemment plus utiles avec l'expérience. Pour l'autre ensemble de constantes donné, le découplage entre les deux composants non multiplié par le$𝜆_ℎ,𝜆_𝑐$ d'ordre $10^-6$ sont très peu variables dans le carré unitaire et très gigantesques à proximité de la perturbation des conditions initiales.
Donc, une solution fermée n'est pas disponible avec les constantes. La question donnée est mal posée et se présente comme une instabilité numérique.
L'ensemble des équations se découpe par $𝜆_ℎ,𝜆_𝑐$.
$(A')$ $\frac{\partial\theta_h}{\partial x}=-\beta_h\theta_h$
$(B')$ $\frac{\partial\theta_c}{\partial x}=-\beta_h\theta_c$
$(C')$->
$(C1)$ $ 𝜆_ℎ\frac{∂^2𝜃_𝑤}{∂𝑥^2}+𝜆_𝑐 𝑉 \frac{∂^2𝜃_𝑤}{∂𝑦^2}=0$
$(C1)$ $−\frac{∂𝜃_h}{∂𝑥}−𝑉\frac{∂𝜃_𝑐}{∂𝑦}=0$
où, $𝛽_ℎ,𝛽_𝑐,𝑉,𝜆_ℎ,𝜆_𝑐$ sont des constantes.
Les conditions aux limites sont:
(JE)
$ \ frac {∂𝜃_𝑤 (0, 𝑦)} {∂𝑥} = \ frac {∂𝜃_𝑤 (1, 𝑦)} {∂𝑥} = \ frac {∂𝜃_𝑤 (𝑥, 0)} {∂𝑦} = \ frac {∂𝜃_𝑤 (𝑥, 1)} {∂𝑦} = 0
Ce sont les conditions aux limites de von Neumann.
Dans Mathematica, il suffit de les saisir de cette manière:
NeumannValue[\[Theta]w[x, y]==0, x == 1 || x == 1 || y == 0 || y == 1];
Cela peut être déduit de la page de message qui est proposée si ceux-ci sont saisis en tant que DirichletConditions
.
Il existe une belle théorie disponible en ligne sur Wolfrom pour estimer les problèmes ou le bien-être du pde: PartialDifferentialEquation .
C'est en quelque sorte un chemin court mais la page de documentation pour le NeumannValue
résout l'équation découplée $ C1 $ avec la pertubation simple disponible. Puisque nous n'avons aucune pertubation. Toutes nos conditions sont nulles sur la frontière. On obtient la solution banale pour $ \ theta_w (x, y) = 0 $ sur le carré entre $ (0,0) $ et $ (1,1) $ .
Mais gardez à l'esprit qu'avec le processus, nous n'obtenons que la solution non homogène. Il y a une solution homogène à ajouter.
Pour présenter la série de Fourier, je me réfère à la page de documentation de DSolve
. De là:
heqn = 0 == D[u[x, t], {x, 2}];
ic = u[x, 0] == 1;
bc = {Derivative[1, 0][u][0, t] == 0,
Derivative[1, 0][u][1, t] == 0};
sol = u[x, t] /. DSolve[{heqn, ic, bc }, u[x, t], {x, t}][[1]]
asol = sol /. {\[Infinity] -> 8} // Activate
Plot3D[asol // Evaluate, {x, 0, 1}, {t, 0, 1}, Exclusions -> None,
PlotRange -> All, AxesLabel -> Automatic]

La solution est DiracDelta[t]
.
Donc rien de vraiment intéressant là-bas. Les conditions aux limites sont remplies. Avec un peu de pertubation, cela ne donnera pas une série de Fourier plus compliquée. DSolve offre quelques exemples. De la série de Fourier, on peut répondre correctement à la première question.
(A ') et (B') sont résolus par des exponentielles qui peuvent être facilement transformées en séries de Fourier.
bh = 0.433; bc = 0.433; \[Lambda]h = 2.33*10^-6; \[Lambda]c =
2.33*10^-6; V = 1;
PDE1 = D[\[Theta]h[x, y], x] + bh*\[Theta]h[x, y] == 0;
PDE2 = D[\[Theta]c[x, y], y] + bc*\[Theta]c[x, y] == 0;
PDE3 = D[\[Theta]h[x, y], x] - V*D[\[Theta]c[x, y], y] == 0;
IC0 = {\[Theta]h[0, y] == 1, \[Theta]c[x, 0] == 0};
(*Random values*)
soli =
NDSolve[{PDE1, PDE2, IC0}, {\[Theta]h, \[Theta]c}, {x, 0, 1}, {y, 0,
1}]

Table[Plot3D[
Evaluate[({\[Theta]h[x, y], \[Theta]c[x, y]} /. soli)[[1, i]]], {x,
0, 1}, {y, 0, 1}, PlotRange -> Full], {i, 1, 2}]

$ \ theta_h (x, y) $ oscille très rapidement sur la frontière et $ \ theta_c (x, y) $ . Par conséquent, toujours dans la solution séparée, il existe une instabilité numérique due à la rigidité de l'accouplement. Seul $ \ theta_c (x, y) $ convient aux conditions initiales mais interfère avec la séparabilité supposée. C'est toujours la double ligne avec pic dans $ \ theta_h (x, y) $ .
Le plus gros problème est la première des conditions initiales.
$$ 𝜃_ℎ (0, 𝑦) = 1, 𝜃_𝑐 (𝑥, 0) = 0 $$
Donc, si pour obtenir une meilleure solution, variez $ 𝜃_ℎ (0, 𝑦) $ ! Faites-le beaucoup plus petit.