Trzy połączone PDE do rozwiązania półanalitycznego / analitycznego

Question

Próbowałem rozwiązać następujące trzy połączone PDE, w których ostatecznym celem jest znalezienie rozkładów $\theta_h, \theta_c$ i $\theta_w$:

$x\in[0,1]$ i $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$$

gdzie, $\beta_h, \beta_c, V, \lambda_h, \lambda_c$są stałymi. Warunki brzegowe to:

$$\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$$

Użytkownik wymiany stosów Matematyki zasugerował mi następujące kroki, które mogą pomóc w rozwiązaniu tego problemu:

1. Przedstaw każdą z trzech funkcji za pomocą dwuwymiarowego szeregu Fouriera
2. Zauważ, że wszystkie równania są liniowe

• Dlatego nie ma sprzężenia częstotliwości
• Tak więc dla każdej pary częstotliwości $\omega_x$, $\omega_y$ będzie rozwiązanie z liniowej kombinacji tylko tych terminów

1. Zastosuj warunki brzegowe bezpośrednio do każdej z trzech serii

• Zauważ, że ze względu na ortogonalność warunek brzegowy musi odnosić się do każdego wyrazu szeregu Fouriera

1. Podłącz szereg Fouriera do PDE i rozwiąż dopasowanie współczynników ( patrz tutaj na przykład w 1D ). Upewnij się, że traktujesz oddzielnie przypadki, w których jedna lub obie częstotliwości są zerowe.
2. Jeśli weźmiesz pod uwagę wszystkie równania dla danej pary częstotliwości, możesz ułożyć je w równanie $M\alpha = 0$, gdzie $\alpha$ są czterema współczynnikami dla tych częstotliwości, i $M$ jest małą rzadką macierzą (np. 12x12), która będzie zależeć tylko od stałych.
3. Dla każdej częstotliwości dozwolone rozwiązania będą znajdować się w zerowej przestrzeni tej macierzy. W przypadku, gdy nie jesteś w stanie analitycznie obliczyć zerowej przestrzeni, nie jest to wielka sprawa - numeryczne obliczenie zerowej przestrzeni jest łatwe, szczególnie w przypadku małych macierzy.

Czy ktoś może mi pomóc w zastosowaniu tych kroków w 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;

Rozwiązanie NDSolve (błędne wyniki)

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

W kierunku rozwiązania możliwego do rozdzielenia

napisałem $\theta_h(x,y) = \beta_h e^{-\beta_h x} \int e^{\beta_h x} \theta_w(x,y) \, \mathrm{d}x$ i $\theta_c(x,y) = \beta_c e^{-\beta_c y} \int e^{\beta_c y} \theta_w(x,y) \, \mathrm{d}y$ i wyeliminowane $\theta_h$ i $\theta_c$z równania (DO). Następnie użyłem ansatz$\theta_w(x,y) = e^{-\beta_h x} f(x) e^{-\beta_c y} g(y)$na tym nowym równaniu. (C), aby go rozdzielić$x$ i $y$składniki. Następnie przy użyciu$F(x) := \int f(x) \, \mathrm{d}x$ i $G(y) := \int g(y) \, \mathrm{d}y$, Otrzymuję następujące dwa równania:

\ 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} z pewną stałą separacji$\mu \in \mathbb{R}$. Nie mogłem jednak kontynuować.

Równanie różniczkowe częściowe całkowitoliczbowe

Eliminowanie $\theta_h, \theta_c$z równania. (C) daje początek całkowaniu równania różniczkowego:

\ 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} \ left (\ lambda_c e ^ {\ beta_c y} \ frac {\ części ^ 2 \ theta_w} {\ części 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}

KOLCE

Dla bc = 4; bh = 2; λc = 0.01; λh = 0.01; V = 2;

Jednak te same parametry, ale V=1ładnie działają.

Materiały referencyjne dla przyszłych użytkowników

Aby zrozumieć ocenę współczynników Fouriera przy użyciu koncepcji minimalizacji najmniejszych kwadratów, którą @bbgodfrey używa w swojej odpowiedzi, przyszli użytkownicy mogą przyjrzeć się tej pracy R. Kelmana (1979). Alternatywnie, ta prezentacja i ten film są również przydatnymi odniesieniami.

Answer 1

Zmiany: Zastąpiono ekspansję na 1 termin ekspansją na okres nokresowy; poprawiona ogólność obliczeń wartości własnych i współczynników; zmieniony i uproszczony kod.

Rozpoczynając od tego zestawu równań, wykonaj następujące czynności, aby uzyskać prawie symboliczne rozwiązanie.

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

Najpierw zamień te równania na równania różniczkowe metodą rozdzielania zmiennych.

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 *)

Mając równania podzielone na równania równań różniczkowych, rozwiąż równania zależne od yz zastosowanymi warunkami brzegowymi. Powstałe wyrażenia, w tym zaangażowane RootSum, są obszerne i dlatego nie są tu przytaczane.

sy = DSolveValue[{eq2y, eq3y, θcy[0] == 0, θwy'[0] == 0}, {θwy[y], θcy[y], θwy'[1]}, 
     {y, 0, 1}] /. C[2] -> 1;

Jest to, oczywiście, problemu wartości własnej, z nietrywialnych rozwiązań tylko dla dyskretnych wartości stałej rozdzielającego sw. Relacja dyspersji dla swjest określona wzorem θwy'[1] == 0. Odpowiednia xzależność jest określana dla każdej wartości własnej przez

sx = DSolveValue[{eq1x, eq3x, θwx'[0] == 0, θwx'[1] == 0, θhx[0] == 1}, 
    {θwx[x], θhx[x]}, {x, 0, 1}];

i to w tym momencie stosuje się niejednorodny warunek brzegowy θhx[0] == 1. Ten wynik jest również zbyt długi, aby go tu odtworzyć.

Następnie numerycznie określ kilka pierwszych (tutaj n = 6) wartości własnych, co wymaga określenia parametrów:

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]

Kilka pierwszych wartości własnych jest szacowanych na podstawie zer wykresu, a następnie obliczanych z dużą dokładnością.

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} *)

i odpowiednie funkcje własne uzyskane przez podłączenie tych wartości swdo sy[1;;2]i 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}]

Po nobliczeniu pierwszych pełnych funkcji własnych określane są ich współczynniki, aby można je było zsumować w celu przybliżenia rozwiązania pierwotnych równań. Odbywa się to metodą najmniejszych kwadratów, ponieważ system ODE nie jest samosprzężony.

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} *)

Jakość dopasowania jest bardzo dobra.

Plot[coef.syn - 1, {y, 0, 1}, AxesLabel -> {y, err}, 
    LabelStyle -> {15, Bold, Black}, ImageSize -> Large]

Na koniec skonstruuj rozwiązanie.

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

Ponieważ to wyprowadzenie jest długie, pokazujemy tutaj, że same równania są spełnione identycznie.

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} *)

Ponadto warunek brzegowy włączony θhjest spełniony w stopniu lepszym niż 0,004%, a warunek brzegowy włączony θcjest spełniony identycznie.

Odpowiednie obliczenia 3D zostały zakończone w 226346 .

Answer 2

Rozwiązanie, które otrzymałem w wersji 12.0.0, wygląda na rzeczywiście niespójne. Porównuję rozwiązanie dość bliskie temu przedstawionemu na stronie dokumentacji NDSolvew sekcji Możliwe problemy -> Równania różniczkowe cząstkowe z przykładem równania Laplace'a z wartościami początkowymi.

Dla podanego układu równań różniczkowych cząstkowych i dla wartości ustawionej tylko jednym mogę użyć NDSolvedo tego wyniku:

Podobieństwo nie polega na rozbieżności, która spada do źródła, ale na rzędzie kolców, które można zobaczyć mniej więcej $x=.3$ i $y=0.3$ dla $𝜃_h$ i $𝜃_c$. To sprzężenie jest jednak naprawdę niefizyczne. Ale jest trochę bardziej pozornie przydatnych informacji dotyczących eksperymentu. Dla drugiego podanego zestawu stałych odsprzężenie między dwoma składowymi nie jest pomnożone przez$𝜆_ℎ,𝜆_𝑐$ zamówienia $10^-6$ bardzo niewiele się różnią w jednostkowej kwadracie i są bardzo gigantyczne w pobliżu zakłóceń z warunków początkowych.

Tak więc rozwiązanie zamknięte nie jest dostępne ze stałymi. Podane pytanie jest źle postawione i przejawia się jako niestabilność liczbowa.

Zbiór równań rozdziela się o $𝜆_ℎ,𝜆_𝑐$.

$(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$

gdzie, $𝛽_ℎ,𝛽_𝑐,𝑉,𝜆_ℎ,𝜆_𝑐$ są stałymi.

Warunki brzegowe to:

(JA)

$ \ frac {∂𝜃_𝑤 (0, 𝑦)} {∂𝑥} = \ frac {∂𝜃_𝑤 (1, 𝑦)} {∂𝑥} = \ frac {∂𝜃_𝑤 (𝑥, 0)} {∂𝑦} = \ frac {∂𝜃_𝑤 (𝑥, 1)} {∂𝑦} = 0

To są warunki brzegowe von Neumanna.

W Mathematica wystarczy wpisać je w ten sposób:

NeumannValue[\[Theta]w[x, y]==0, x == 1 || x == 1 || y == 0 || y == 1];

Można to wywnioskować ze strony wiadomości, która jest oferowana, jeśli zostaną wprowadzone jako DirichletConditions.

Wolfrom udostępnia w Internecie fajną teorię do szacowania problemów lub dobrego zachowania pde: PartialDifferentialEquation .

Jest to w pewnym sensie krótka trasa, ale strona z dokumentacją NeumannValuerozwiązuje rozdzielone równanie $ C1 $ z dostępną prostą pertubacją. Ponieważ nie mamy pertubacji. Wszystkie nasze warunki są równe zeru na granicy. Otrzymujemy banalne rozwiązanie dla $ \ theta_w (x, y) = 0 $ na kwadracie między $ (0,0) $ a $ (1,1) $ .

Pamiętaj jednak, że w tym procesie otrzymujemy tylko niejednorodne rozwiązanie. Należy dodać jednorodny roztwór.

Aby przedstawić szereg Fouriera, odsyłam do strony dokumentacji DSolve. Stamtąd:

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]

Rozwiązaniem jest DiracDelta[t].

Więc nie ma tam nic naprawdę interesującego. Warunki brzegowe są spełnione. Przy pewnej pertubacji da to bardziej skomplikowany szereg Fouriera. DSolve oferuje kilka przykładów. Z szeregu Fouriera na pierwsze pytanie można odpowiedzieć poprawnie.

(A ') i (B') są rozwiązywane przez wykładniki, które można wygodnie przekształcić w szeregi Fouriera.

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) $ oscyluje bardzo szybko na granicy, a $ \ theta_c (x, y) $ . Dlatego nadal w rozdzielonym rozwiązaniu występuje niestabilność numeryczna spowodowana sztywnością sprzęgła. Tylko $ \ theta_c (x, y) $ pasuje do warunków początkowych, ale koliduje z założoną rozdzielnością. Nadal jest to podwójny wiersz ze skokiem w $ \ theta_h (x, y) $ .

Największym problemem jest pierwszy z warunków początkowych.

$$ 𝜃_ℎ (0, 𝑦) = 1, 𝜃_𝑐 (𝑥, 0) = 0 $$

Więc jeśli chcesz uzyskać ładniejsze rozwiązanie, zmień $ 𝜃_ℎ (0, 𝑦) $ ! Spraw, aby był znacznie mniejszy.