반 분석 / 분석적으로 풀어야 할 3 개의 결합 된 PDE

Nov 30 2020

나는 최종 목표가 분포를 찾는 다음 세 가지 PDE를 해결하려고 노력해 왔습니다. $\theta_h, \theta_c$$\theta_w$:

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

어디, $\beta_h, \beta_c, V, \lambda_h, \lambda_c$상수입니다. 경계 조건은 다음과 같습니다.

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

Mathematics 스택 교환의 한 사용자가이 문제를 해결하는 데 도움이 될 수있는 다음 단계를 제안했습니다.

  1. 2D 푸리에 급수를 사용하여 세 가지 함수를 각각 나타냅니다.
  2. 모든 방정식이 선형인지 관찰하십시오.
  • 따라서 주파수 결합이 없습니다.
  • 따라서 모든 주파수 쌍에 대해 $\omega_x$, $\omega_y$ 해당 항만 선형 조합에서 솔루션이 있습니다.
  1. 세 시리즈 각각에 직접 경계 조건 적용
  • 직교성에 따라 경계 조건은 푸리에 급수의 각 항에 적용되어야합니다.
  1. 푸리에 급수를 PDE에 연결하고 계수 일치를 풉니 다 ( 예 : 1D 참조 ). 주파수 중 하나 또는 둘 다가 0 인 경우를 별도로 처리해야합니다.
  2. 주어진 주파수 쌍에 대한 모든 방정식을 고려하면 방정식으로 배열 할 수 있습니다. $M\alpha = 0$, 어디 $\alpha$ 그 주파수에 대한 푸리에 계수이고, $M$ 상수에만 의존하는 작은 희소 행렬 (12x12와 같은)입니다.
  3. 각 주파수에 대해 허용되는 솔루션은 해당 행렬의 Null 공간에 있습니다. 널 공간을 분석적으로 해결할 수없는 경우에는 큰 문제가 아닙니다. 특히 작은 행렬의 경우 널 공간을 수치 적으로 계산하는 것이 쉽습니다.

누군가 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;

NDSolve 솔루션 (잘못된 결과)

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

분리 가능한 솔루션을 향하여

나는 썼다 $\theta_h(x,y) = \beta_h e^{-\beta_h x} \int e^{\beta_h x} \theta_w(x,y) \, \mathrm{d}x$$\theta_c(x,y) = \beta_c e^{-\beta_c y} \int e^{\beta_c y} \theta_w(x,y) \, \mathrm{d}y$ 그리고 제거 $\theta_h$$\theta_c$식에서. (씨). 그런 다음 ansatz를 사용했습니다.$\theta_w(x,y) = e^{-\beta_h x} f(x) e^{-\beta_c y} g(y)$이 새로운 Eq. (C) 분리$x$$y$구성 요소. 그런 다음 사용$F(x) := \int f(x) \, \mathrm{d}x$$G(y) := \int g(y) \, \mathrm{d}y$, 다음 두 가지 방정식을 얻습니다.

\ 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 '' '-2V \ 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} ( 일부 분리 상수 포함)$\mu \in \mathbb{R}$. 그러나 나는 더 이상 진행할 수 없었다.

부분 정수 미분 방정식

제거 $\theta_h, \theta_c$Eq. (C) 부분 적분 미분 방정식을 생성합니다.

\ 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} \ 왼쪽 (\ 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}

스파이크

에 대한 bc = 4; bh = 2; λc = 0.01; λh = 0.01; V = 2;

그러나 동일한 매개 변수이지만 V=1잘 작동합니다.

향후 사용자를위한 참고 자료

@bbgodfrey 가 답변에서 사용 하는 최소 제곱의 최소화 개념을 사용하여 푸리에 계수의 평가를 이해하기 위해 향후 사용자는 R. Kelman (1979) 의이 논문을 볼 수 있습니다 . 또는 이 프레젠테이션 과이 비디오도 유용한 참고 자료입니다.

답변

3 bbgodfrey Dec 06 2020 at 09:14

편집 : 1 항 확장을 n 항 확장으로 대체했습니다. 고유 값 및 계수 계산의 일반성 개선; 재정렬되고 단순화 된 코드.

이 방정식 세트로 시작하여 거의 상징적 인 솔루션을 얻으려면 다음과 같이 진행하십시오.

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

먼저 변수 분리 방법을 사용하여 이러한 방정식을 ODE로 변환합니다.

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

방정식을 ODE로 분리하여 경계 조건이 적용된 y 종속 방정식을 풉니 다. 를 포함하는 결과 표현식 RootSum은 길기 때문에 여기에서 재현되지 않습니다.

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

물론 이것은 분리 상수의 이산 값에 대해서만 사소하지 않은 솔루션이있는 고유 값 문제입니다 sw. 에 대한 분산 관계는 sw로 제공됩니다 θwy'[1] == 0. 해당 x종속성은 다음과 같이 각 고유 값에 대해 결정됩니다.

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

이 시점에서 비균질 경계 조건 θhx[0] == 1,가 적용됩니다. 이 결과도 여기에서 재현하기에는 너무 길다.

다음으로, n = 6매개 변수를 지정해야하는 처음 몇 개의 고유 값 (여기서는 )을 수치 적으로 결정합니다 .

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]

처음 몇 개의 고유 값은 플롯의 0에서 추정 된 다음 높은 정확도로 계산됩니다.

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

및이 값 swsy[1;;2]및에 대입하여 얻은 해당 고유 함수 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}]

첫 번째 n완전한 고유 함수가 계산되면 다음으로 계수가 결정되어 원래 방정식에 대한 해를 근사하기 위해 합산 될 수 있습니다. 이것은 ODE 시스템이 자기 결합이 아니기 때문에 최소 제곱 법에 의해 수행됩니다.

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

핏의 품질이 매우 좋습니다.

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

마지막으로 솔루션을 구성하십시오.

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

이 도출은 길기 때문에 방정식 자체가 동일하게 충족된다는 것을 여기서 보여줍니다.

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

또한 경계 조건 on θh은 0.004 % 이상으로 만족하며, 경계 조건 θc은 동일하게 만족합니다.

해당 3D 계산은 226346 에서 완료되었습니다 .

2 SteffenJaeschke Dec 06 2020 at 00:34

버전 12.0.0에서 얻은 솔루션은 실제로 일관성이 없어 보입니다. NDSolve가능한 문제-> 편미분 방정식 섹션 의 문서 페이지에 표시된 솔루션 과 초기 값이있는 라플라스 방정식의 예를 비교합니다.

주어진 편미분 방정식 시스템과 하나만으로 설정된 값 NDSolve에 대해이 결과에 사용할 수 있습니다.

유사성은 원점으로 떨어지는 발산이 아니라 약에서 볼 수있는 스파이크 열입니다. $x=.3$$y=0.3$ ...에 대한 $𝜃_h$$𝜃_c$. 이 결합은 실제로 비 물리적입니다. 그러나 실험에 대해 좀 더 유용한 정보가 있습니다. 주어진 다른 상수 세트의 경우 두 구성 요소 사이의 디커플링은$𝜆_ℎ,𝜆_𝑐$ 주문 $10^-6$ 단위 광장에서 거의 변하지 않으며 초기 조건의 교란에 매우 가깝습니다.

따라서 닫힌 솔루션은 상수와 함께 사용할 수 없습니다. 주어진 질문은 잘못된 자세로 수치 적 불안정으로 나타납니다.

방정식 세트는 다음과 같이 분리됩니다. $𝜆_ℎ,𝜆_𝑐$.

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

어디, $𝛽_ℎ,𝛽_𝑐,𝑉,𝜆_ℎ,𝜆_𝑐$ 상수입니다.

경계 조건은 다음과 같습니다.

(나는)

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

이것이 폰 노이만 경계 조건입니다.

Mathematica에서는 다음과 같이 입력하면 충분합니다.

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

로 입력 된 경우 제공되는 메시지 페이지에서 유추 할 수 있습니다 DirichletConditions.

pde : PartialDifferentialEquation 의 문제 또는 웰빙을 평가하기 위해 Wolfrom에서 온라인으로 사용할 수있는 몇 가지 좋은 이론이 있습니다 .

어떻게 든 짧은 경로이지만 설명서 페이지 는 몇 가지 간단한 섭동을 사용 NeumannValue하여 분리 방정식 $ C1 $ 을 해결합니다 . 우리는 삽관이 없기 때문에. 우리의 모든 조건은 경계에서 0입니다. $ (0,0) $$ (1,1) $ 사이의 정사각형 에서 $ \ theta_w (x, y) = 0 $ 에 대한 평범한 해를 얻습니다 .

그러나 우리는 비균질적인 해결책만을 얻는 과정을 명심하십시오. 추가 할 균질 솔루션이 있습니다.

푸리에 시리즈를 소개하기 위해 문서 페이지를 참조합니다 DSolve. 거기에서:

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]

해결책은 DiracDelta[t].

그래서 거기에는 정말 흥미로운 것이 없습니다. 경계 조건이 충족됩니다. 약간의 삽관으로이 woult는 더 복잡한 푸리에 급수를 제공합니다. DSolve 는 몇 가지 예를 제공합니다. 푸리에 시리즈에서 첫 번째 질문에 대한 올바른 답을 찾을 수 있습니다.

(A ')와 (B')는 편안하게 푸리에 급수로 변환 될 수있는 지수로 풀립니다.

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) $ 는 경계와 $ \ theta_c (x, y) $ 에서 매우 빠르게 진동 합니다. 따라서 여전히 분리 된 솔루션에는 커플 링의 강성으로 인한 수치 적 불안정성이 있습니다. 단지 $ \ theta_c (X, Y) $ 적합한 가정 된 분리도 초기 조건이지만 방해한다. $ \ theta_h (x, y) $에 스파이크가있는 이중 행입니다 .

가장 큰 문제는 첫 번째 초기 조건입니다.

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

따라서 더 좋은 솔루션을 얻으려면 $ 𝜃_ℎ (0, 𝑦) $ ! 훨씬 작게 만드십시오.