Drei gekoppelte PDEs müssen semi-analytisch / analytisch gelöst werden
Ich habe versucht, die folgenden drei gekoppelten PDEs zu lösen, bei denen das endgültige Ziel darin besteht, die Verteilungen zu finden $\theta_h, \theta_c$ und $\theta_w$::
$x\in[0,1]$ und $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$$
wo, $\beta_h, \beta_c, V, \lambda_h, \lambda_c$sind Konstanten. Die Randbedingungen sind:
$$\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$$
Ein Benutzer von Mathematics Stack Exchange schlug mir die folgenden Schritte vor, um dieses Problem zu lösen:
- Stellen Sie jede der drei Funktionen mithilfe von 2D-Fourier-Reihen dar
- Beachten Sie, dass alle Gleichungen linear sind
- Somit gibt es keine Frequenzkopplung
- Also für jedes Frequenzpaar $\omega_x$, $\omega_y$ Es wird eine Lösung aus einer linearen Kombination nur dieser Begriffe geben
- Wenden Sie die Randbedingungen direkt auf jede der drei Serien an
- Beachten Sie, dass durch Orthogonalität die Randbedingung für jeden Term der Fourier-Reihe gelten muss
- Stecken Sie die Fourier-Reihe in die PDE und lösen Sie die Koeffizientenanpassung ( siehe hier zum Beispiel in 1D ). Stellen Sie sicher, dass Sie die Fälle, in denen eine oder beide Frequenzen Null sind, separat behandeln.
- Wenn Sie alle Gleichungen für ein bestimmtes Frequenzpaar berücksichtigen, können Sie sie in einer Gleichung anordnen $M\alpha = 0$, wo $\alpha$ sind Fourier-Koeffizienten für diese Frequenzen und $M$ ist eine kleine, spärliche Matrix (etw wie 12x12), die nur von den Konstanten abhängt.
- Für jede Frequenz befinden sich die zulässigen Lösungen im Nullraum dieser Matrix. Falls Sie nicht in der Lage sind, den Nullraum analytisch zu lösen, ist dies keine große Sache - die numerische Berechnung des Nullraums ist insbesondere für kleine Matrizen einfach.
Kann mir jemand bei der Anwendung dieser Schritte in Mathematica helfen?
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;
NDSolve-Lösung (falsche Ergebnisse)
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}]
Auf dem Weg zu einer trennbaren Lösung
Ich schrieb $\theta_h(x,y) = \beta_h e^{-\beta_h x} \int e^{\beta_h x} \theta_w(x,y) \, \mathrm{d}x$ und $\theta_c(x,y) = \beta_c e^{-\beta_c y} \int e^{\beta_c y} \theta_w(x,y) \, \mathrm{d}y$ und beseitigt $\theta_h$ und $\theta_c$aus Gl. (C). Dann habe ich den Ansatz benutzt$\theta_w(x,y) = e^{-\beta_h x} f(x) e^{-\beta_c y} g(y)$auf dieser neuen Gl. (C) um es zu trennen$x$ und $y$Komponenten. Dann weiter$F(x) := \int f(x) \, \mathrm{d}x$ und $G(y) := \int g(y) \, \mathrm{d}y$Ich erhalte die folgenden zwei Gleichungen:
\ 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' '+ \ left ((\ lambda_c \ beta_c - 1) V \ beta_c + \ mu \ right) G' + V \ beta_c ^ 2 G & = & 0, \ end {eqnarray} mit einer gewissen Trennungskonstante$\mu \in \mathbb{R}$. Ich konnte jedoch nicht weiter vorgehen.
Eine partielle Integro-Differentialgleichung
Beseitigen $\theta_h, \theta_c$aus der Gl. (C) führt zu einer partio-integralen Differentialgleichung:
\ begin {eqnarray} 0 & = & e ^ {- \ beta_h x} \ left (\ lambda_h e ^ {\ beta_h x} \ frac {\ partielle ^ 2 \ theta_w} {\ partielle 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} \ links (\ lambda_c e ^ {\ beta_c y} \ frac {\ partiell ^ 2 \ theta_w} {\ partiell 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}
SPIKES
Zum bc = 4; bh = 2; λc = 0.01; λh = 0.01; V = 2;
Allerdings V=1funktionieren die gleichen Parameter aber gut.
Einige Referenzmaterialien für zukünftige Benutzer
Um die Bewertung von Fourier-Koeffizienten unter Verwendung des Konzepts der Minimierung kleinster Quadrate zu verstehen, das @bbgodfrey in seiner Antwort verwendet, können zukünftige Benutzer dieses Papier von R. Kelman (1979) betrachten. Alternativ sind diese Präsentation und dieses Video auch nützliche Referenzen.
Antworten
Änderungen: 1-Term-Erweiterung durch n-Term-Erweiterung ersetzt; verbesserte Allgemeinheit der Eigenwert- und Koeffizientenberechnungen; neu geordneter und vereinfachter Code.
Beginnen Sie mit diesem Satz von Gleichungen wie folgt, um eine fast symbolische Lösung zu erhalten.
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])
Konvertieren Sie diese Gleichungen zunächst durch die Methode der Trennung von Variablen in ODEs.
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 *)
Lösen Sie mit den in ODEs getrennten Gleichungen die y-abhängigen Gleichungen mit den angewendeten Randbedingungen. Die daraus resultierenden Ausdrücke RootSumsind langwierig und werden daher hier nicht wiedergegeben.
sy = DSolveValue[{eq2y, eq3y, θcy[0] == 0, θwy'[0] == 0}, {θwy[y], θcy[y], θwy'[1]},
{y, 0, 1}] /. C[2] -> 1;
Dies ist natürlich ein Eigenwertproblem mit nichttrivialen Lösungen nur für diskrete Werte der Trennungskonstante sw. Die Dispersionsrelation für swist gegeben durch θwy'[1] == 0. Die entsprechende xAbhängigkeit wird für jeden Eigenwert durch bestimmt
sx = DSolveValue[{eq1x, eq3x, θwx'[0] == 0, θwx'[1] == 0, θhx[0] == 1},
{θwx[x], θhx[x]}, {x, 0, 1}];
und an diesem Punkt wird die inhomogene Randbedingung θhx[0] == 1angewendet. Dieses Ergebnis ist auch zu lang, um es hier wiederzugeben.
Bestimmen Sie als Nächstes numerisch die ersten mehreren (hier n = 6) Eigenwerte, für die die Parameter angegeben werden müssen:
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]
Die ersten mehreren Eigenwerte werden aus den Nullen des Diagramms geschätzt und dann mit hoher Genauigkeit berechnet.
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} *)
und die entsprechenden Eigenfunktionen, die durch Einstecken dieser Werte swin sy[1;;2]und erhalten werden 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}]
Wenn die ersten nvollständigen Eigenfunktionen berechnet sind, werden als nächstes ihre Koeffizienten bestimmt, so dass sie summiert werden können, um die Lösung an die ursprünglichen Gleichungen anzunähern. Dies geschieht durch kleinste Quadrate, da das ODE-System nicht selbstadjunkt ist.
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} *)
Die Qualität der Passform ist sehr gut.
Plot[coef.syn - 1, {y, 0, 1}, AxesLabel -> {y, err},
LabelStyle -> {15, Bold, Black}, ImageSize -> Large]
Schließlich konstruieren Sie die Lösung.
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}]
Da diese Ableitung langwierig ist, zeigen wir hier, dass die Gleichungen selbst identisch erfüllt sind.
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} *)
Darüber hinaus ist die Randbedingung ein θhbesser als 0,004% erfüllt, und die Randbedingung ein θcist identisch erfüllt.
Die entsprechende 3D-Berechnung wurde um 226346 abgeschlossen .
Die Lösung, die ich mit Version 12.0.0 bekomme, sieht in der Tat inkonsistent aus. Ich vergleiche die Lösung ziemlich nahe an der auf der Dokumentationsseite für NDSolveim Abschnitt Mögliche Probleme -> Partielle Differentialgleichungen gezeigten mit dem Beispiel für die Laplace-Gleichung mit Anfangswerten.
Für das angegebene partielle Differentialgleichungssystem und für den nur mit einem eingestellten Wert kann ich NDSolvefür dieses Ergebnis Folgendes verwenden:
Die Ähnlichkeit ist nicht die Divergenz, die zum Ursprung abfällt, sondern die Reihe von Spitzen, die ungefähr zu sehen sind $x=.3$ und $y=0.3$ zum $𝜃_h$ und $𝜃_c$. Diese Kopplung ist jedoch wirklich unphysisch. Das Experiment enthält jedoch einige weitere scheinbar nützliche Informationen. Für den anderen gegebenen Satz von Konstanten wird die Entkopplung zwischen den beiden Komponenten nicht mit der multipliziert$𝜆_ℎ,𝜆_𝑐$ der Ordnung $10^-6$ sind sehr wenig variabel im Einheitsquadrat und sehr gigantisch in der Nähe der Störung von den Anfangsbedingungen.
Daher ist mit den Konstanten keine geschlossene Lösung verfügbar. Die gegebene Frage ist schlecht gestellt und zeigt sich als numerische Instabilität.
Der Gleichungssatz entkoppelt sich um $𝜆_ℎ,𝜆_𝑐$.
$(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$
wo, $𝛽_ℎ,𝛽_𝑐,𝑉,𝜆_ℎ,𝜆_𝑐$ sind Konstanten.
Die Randbedingungen sind:
(ICH)
$ \ frac {∂𝜃_𝑤 (0, 𝑦)} {∂𝑥} = \ frac {∂𝜃_𝑤 (1, 𝑦)} {∂𝑥} = \ frac {∂𝜃_𝑤 (𝑥, 0)} {∂𝑦} = \ frac {∂𝜃_𝑤 (𝑥, 1)} {∂𝑦} = 0
Dies sind von Neumann-Randbedingungen.
In Mathematica reicht es aus, sie folgendermaßen einzugeben:
NeumannValue[\[Theta]w[x, y]==0, x == 1 || x == 1 || y == 0 || y == 1];
Dies kann aus der angebotenen Nachrichtenseite abgeleitet werden, wenn diese als eingegeben werden DirichletConditions.
Es gibt eine nette Theorie online von Wolfrom, um die Probleme oder das Verhalten der PDE abzuschätzen: PartialDifferentialEquation .
Es ist irgendwie ein kurzer Weg, aber die Dokumentationsseite für NeumannValuelöst die entkoppelte Gleichung $ C1 $ mit einer einfachen Pertubation. Da haben wir keine Störung. Alle unsere Bedingungen sind an der Grenze Null. Wir erhalten die banale Lösung für $ \ theta_w (x, y) = 0 $ auf dem Quadrat zwischen $ (0,0) $ und $ (1,1) $ .
Denken Sie jedoch daran, dass wir nur die inhomogene Lösung erhalten. Es muss eine homogene Lösung hinzugefügt werden.
Zur Einführung der Fourier-Reihe verweise ich auf die Dokumentationsseite von DSolve. Von dort:
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]
Die Lösung ist DiracDelta[t].
Also nichts wirklich interessantes da. Die Randbedingungen sind erfüllt. Mit etwas Pertubation ergibt dieses Ergebnis eine kompliziertere Fourier-Reihe. DSolve bietet einige Beispiele. Aus der Fourier-Reihe kann die erste Frage richtig beantwortet werden.
(A ') und (B') werden durch Exponentiale gelöst, die bequem in Fourier-Reihen umgewandelt werden können.
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) $ schwingt sehr schnell an der Grenze und $ \ theta_c (x, y) $ . Daher besteht immer noch in der getrennten Lösung eine numerische Instabilität aufgrund der Steifheit der Kupplung. Nur $ \ theta_c (x, y) $ entspricht den Anfangsbedingungen, beeinträchtigt jedoch die angenommene Trennbarkeit. Es ist immer noch die doppelte Reihe mit Spitze in $ \ theta_h (x, y) $ .
Das größte Problem ist die erste der Anfangsbedingungen.
$$ 𝜃_ℎ (0, 𝑦) = 1, 𝜃_𝑐 (𝑥, 0) = 0 $$
Wenn Sie also eine schönere Lösung erhalten möchten, variieren Sie $ 𝜃_ℎ (0, 𝑦) $ ! Mach es viel kleiner.