# Solving a PDE (2D Laplacian) coupled with an ODE

I have asked this question before, but this is my new attempt and so instead of cluttering the previous one, I am making a new post. I am trying to analytically solve a PDE ($$\nabla^2 T(x,y)=0$$) coupled with an ODE. The PDE is subjected to the following boundary conditions:

$$\frac{\partial T(0,y)}{\partial x}=\frac{\partial T(L,y)}{\partial x}=0 \tag 1$$

$$\frac{\partial T(x,0)}{\partial y}=\gamma \tag 2$$

$$\frac{\partial T(x,l)}{\partial y}=\beta (T(x,l)-t) \tag 3$$

where $$t$$ is governed by the ODE:

$$\frac{\partial t}{\partial x}+\alpha(t-T(x,l))=0 \tag 4$$

subjected to $$t(x=0)=0$$. I am trying separation of variables. I manipulated $$(4)$$ to express $$t$$ as $$t=\alpha e^{-\alpha x}\Bigg(\int_0^x e^{\alpha s }T(s,l)\mathrm{d}s\Bigg)$$ and substituted in $$(3)$$ while applying the 3rd b.c.

My attempt is (I must acknowledge Bill Watts here as I have used methods I learned from his answer on MMA SE):

pde = D[T[x, y], x, x] + D[T[x, y], y, y] == 0

(*product form*)
T[x_, y_] = X[x] Y[y]
pde/T[x, y] // Expand
xeq = X''[x]/X[x] == -a^2
DSolve[xeq, X[x], x] // Flatten
X[x_] = X[x] /. % /. {C[1] -> c1, C[2] -> c2}
yeq = Y''[y]/Y[y] == a^2
DSolve[yeq, Y[y], y] // Flatten
Y[y_] = (Y[y] /. % /. {C[1] -> c3, C[2] -> c4})

T[x_, y_] = Xp[x] + Yp[y]
xpeq = Xp''[x] == b
DSolve[xpeq, Xp[x], x] // Flatten
Xp[x_] = Xp[x] /. % /. {C[1] -> c5, C[2] -> c6}
ypeq = Yp''[y] + b == 0
DSolve[ypeq, Yp[y], y] // Flatten
Yp[y_] = Yp[y] /. % /. {C[1] -> 0, C[2] -> c7}

T[x_, y_] = X[x] Y[y] + Xp[x] + Yp[y]
pde // FullSimplify

(*Applying the first and second b.c.*)
(D[T[x, y], x] /. x -> 0) == 0
c6 = 0
c2 = 0
c1 = 1

(D[T[x, y], x] /. x -> L) == 0
b = 0
a = (n π)/L
$Assumptions = n ∈ Integers (*Applying the third b.c.*) (D[T[x, y], y] /. y -> 0) == γ c4 = c4 /. Solve[Coefficient[%[[1]], Cos[(π n x)/L]] == 0, c4][[1]] c7 = c7 /. Solve[c7 == γ, c7][[1]] T[x, y] // Collect[#, c3] & (*now splitting T[x,y] into two parts*) T[x, y] /. n -> 0 T0[x_, y_] = 2 c3 + c5 + y γ /. c5 -> 0 Tn[x_, y_] = T[x, y] - T0[x, y] // Simplify (*applying the fourth b.c. to each part individually and using orthogonality*) bcfn0 = (D[T0[x, y], y] /. y -> l) == β (T0[x, l] - α E^(-α x) Integrate[E^(α s) T0[s, l], {s, 0, x}]) Integrate[bcfn0[[1]], {x, 0, L}] == Integrate[bcfn0[[2]], {x, 0, L}] Solve[%, c3] c3 = c3 /. %[[1]] bcfn = (D[Tn[x, y], y] /. y -> l) == β (Tn[x, l] - α E^(-α x) Integrate[E^(α s) Tn[s, l], {s, 0, x}]) Solve[Integrate[bcfn[[1]]*Cos[(n*Pi*x)/L], {x, 0, L}] == Integrate[bcfn[[2]]*Cos[(n*Pi*x)/L], {x, 0, L}], c5]; c5 = c5 /. %[[1]];//FullSimplify T0[x_, y_] = T0[x, y] // Simplify Tn[x_, y_] = Tn[x, y] // Simplify  Now we declare some constants and compile the functions α = 62.9/2; β = 1807/390; γ = 3091.67/390; L = 0.060; l = 0.003; T[x_, y_, mm_] := T0[x, y] + Sum[Tn[x, y], {n, 1, mm}] Plot[{Evaluate[T[x, 0, 10]], Evaluate[T[x, l/2, 10]], Evaluate[T[x, l, 10]]}, {x, 0, L}]  The plot results are extremely ambiguous. The solution is not even converging (as I increase the number of terms , the T value keeps increasing). I cannot figure out what I have done wrong. Since the $$T$$ results are completely out, I have not calculated $$t$$. I cannot figure out what I have done wrong. • Are you sure that the t[0] == 0 condition is correct? – bbgodfrey Aug 27 '20 at 2:28 • If it helps, then this condition can be taken as$t=t_i$Then the$t$differential equation evaluates to $$t=\alpha e^{-\alpha x}\Bigg(\int_0^x e^{\alpha s }T(s,l)\mathrm{d}s+\frac{t_i}{\alpha}\Bigg)$$$t_i\$ is a constant whose value can be taken as 300 – Indrasis Mitra Aug 27 '20 at 3:11

I can fix the problem of your solution increasing with increasing n, but that will not give you a solution. Rather than copy your entire solution I will start where I think the problem starts.

You have

T0[x_, y_] = 2 c3 + c5 + y γ /. c5 -> 0


Change that to

T0[x_, y_] = 2 c3 + c5 + y γ /. c3 -> 0
(*c5 + γ y*)


Then

Tn[x_, y_] = T[x, y] - T0[x, y] // FullSimplify
(*2 c3 Cos[(π n x)/L] Cosh[(π n y)/L]*)


In your case you were carrying an extra constant term c5 with Tn which was being added for every term in your sum which is why your solution increased with each term. In my case I carry c5 as the constant term, but only with T0. The changes below will require changing solving for c5 with bcf0 and solving for c3 with bcfn.

This next problem I fear is insurmountable with the calculation of bcfn0.

bcfn0 = (D[T0[x, y], y] /. y -> l) == β (T0[x, l] - α E^(-α x) Integrate[E^(α s) T0[s, l], {s, 0, x}]) // FullSimplify
(*γ E^(α x) == β (c5 + γ l)*)


Examining this result, it is obvious that there is no constant value c5 can take to satisfy this equation.

Furthermore, with the new Tn the orthogonality equation will result in c3 = 0. This means that T will have no x dependence, which when you think about it makes sense, if T is to satisfy Laplace eq and have x derivatives equal to zero at both ends in the x direction.

If T has no x dependence, then its derivatives can also have no x dependence, but with the y derivative of T depending on t which has x dependence, we have a problem.

• Hi Bill, thanks. I had actually tried with c3->0, to reach a Tn=0 scenario. But I had no idea why in my present solution I had no convergence. But this became clear from your explanation in this answer. I think the problem lies with the convective type b.c. which seems physically correct, but is causing problems. – Indrasis Mitra Aug 28 '20 at 4:16
• Also, have you had any success with the cross-flow problem ? – Indrasis Mitra Aug 28 '20 at 4:19
• I had the exact same problem with the 3D problem. Tn = 0 and T0 not satisfying the bc. – Bill Watts Aug 28 '20 at 4:24
• @IndrasisMitra I do not believe that this problem has a solution. You might consider posting it to Mathematics.StackExchange, because it looks more like a math than a Mathematica problem. – bbgodfrey Aug 28 '20 at 18:34
• I will look at it when I get the chance to to do so in detail. I think this paper does a better job on details, than the other one, but we'll see. As far as the other paper, if you don't have the Stehfest algorithm, you can find it at library.wolfram.com/infocenter/MathSource/2691 if you haven't discovered it already. – Bill Watts Sep 1 '20 at 0:37