# Laplace's equation with mixed boundary condition using separation of variables

The equation and boundary condition are defined in the picture where $$T_1$$, $$T_2$$, $$T_3$$, $$k$$ and $$h$$ are constant value.

I am trying to use variable separation to solve the problem. If $$T(x,y)$$ could be separated to $$X(x)Y(y)$$, then I can get the general solution of the pde:

$$X=A\cos(\lambda x)+B\sin(\lambda x)$$

$$Y=Ce^{- \lambda y}+De^{ \lambda y}$$

$$T(x,y)=(A\cos( \lambda x)+B\sin(\lambda x))(Ce^{- \lambda y}+De^{\lambda y})$$

$$\lambda$$ is the separation constant.

Is there anyway that I can use Mathematica to solve for the $$A,B,C,D$$ and $$\lambda$$?

Thank you for helping me!

• What exactly are you trying to do? Implement separation of variable in Mathematica, or just solve the problem symbolically with Mathematica? Jul 22, 2020 at 7:24
• DSolve will use separation of variables internally for this pde. You do not have to implement it yourself. Jul 22, 2020 at 9:37
• The primary goal is to solve the pde symbolically, Jul 22, 2020 at 15:27
• A couple of points: (1) By redefining $T \to T - T_2$, you can make two of the BCs homogeneous. (2) It's not clear to me that you shouldn't also be including terms of the form $T(x,y)=(A\cos( \lambda y)+B\sin(\lambda y))(Ce^{- \lambda x}+De^{\lambda x})$ as well. Jul 22, 2020 at 16:05

Sorry, I found minor issue actually after I verified with Maple. When breaking the problem into 4 problems, to make life easier for DSolve, the top BC, should be k*Derivative[0, 1][T][x, b] == h*(T[x, b] - T3) when non homogeneous, but for the homogeneous case, it should be k*Derivative[0, 1][T][x, b] == h*(T[x, b]) and not k*Derivative[0, 1][T][x, b] == 0 as I had it below. i.e. only T3 should be set to zero.

I noticed this when I compared Maple's solution to Mathematica's and found very small difference in numerical values.

But now DSolve now can't solve it any more :(

So will post Maple's solution below, and post the corrected Mathematica solution, which breaks the problems into 4, but now it does not solve it. So may be you can try NDSolve in Mathematica for this. May be in V 12.2 DSolve can do it.

## Maple solution

restart;
pde := diff(T(x, y),x$$2)+diff(T(x,y),y$$2)=0;
bc := T(0, y) = T1, T(a, y) = T2, T(x, 0) =T2, k*D[2](T)(x,b)= h*(T(x, b) - T3);
sol1:=simplify(pdsolve([pde, bc], T(x, y)) assuming a>0,b>0);


which gives

T(x,y) = (Sum(-2*(-(-T1+T2)*(-Pi*k*n+a*h)*exp(n*Pi*(2*b-y)/a)+a*((T2-T3)*(-1)^n
-T1+T3)*h*exp(n*Pi*(b-y)/a)-a*((T2-T3)*(-1)^n-T1+T3)*h*exp(n*Pi*(y+b)/a)+exp(n/
a*Pi*y)*(-T1+T2)*(Pi*k*n+a*h))*sin(n/a*Pi*x)/Pi/n/((-Pi*k*n+a*h)*exp(2*n/a*Pi*b
)-Pi*n*k-a*h),n = 1 .. infinity)*a+T1*a+x*(-T1+T2))/a


## Corrected Mathematica solution

ClearAll[T, x, y, T1, T2, T3, k, h];
pde = Laplacian[T[x, y], {x, y}] == 0;
bc = {T[0, y] == 0, T[a, y] == 0, T[x, 0] == 0, k*Derivative[0, 1][T][x, b] - h*T[x, b] == -h*T3}
sol1 = DSolve[{pde, bc}, T[x, y], {x, y}, Assumptions -> {a > 0, b > 0}]


bc = {T[0, y] == 0, T[a, y] == 0, T[x, 0] == T2, k*Derivative[0, 1][T][x, b] - h*T[x, b] == 0}
sol2 = DSolve[{pde, bc}, T[x, y], {x, y},   Assumptions -> {a > 0, b > 0}]


bc = {T[0, y] == 0, T[a, y] == T2, T[x, 0] == 0,   k*Derivative[0, 1][T][x, b] - h*T[x, b] == 0}
sol3 = DSolve[{pde, bc}, T[x, y], {x, y},   Assumptions -> {a > 0, b > 0}]

(* no solution *)

bc = {T[0, y] == T1, T[a, y] == 0, T[x, 0] == 0, k*Derivative[0, 1][T][x, b] - h*T[x, b] == 0}
sol4 = DSolve[{pde, bc}, T[x, y], {x, y},   Assumptions -> {a > 0, b > 0}]

(* no solution *)


So can't add them, since 2 solutions could not be found. need all 4.

## Original answer below (but contain error in BC)

You can solve it as follows.

DSolve can not handle more than one edge in Laplacian in 2D being non homogeneous at same time.

Due to linearity, the problem is broken into 4 problems, were one edge is non homogeneous at time. Then the 4 solutions are added.

ClearAll[T, x, y, T1, T2, T3, k, h];
pde = Laplacian[T[x, y], {x, y}] == 0;
bc = {T[0, y] == T1, T[a, y] == T2, T[x, 0] == T2, k*Derivative[0, 1][T][x, b] == h*(T[x, b] - T3)}
sol1 = DSolve[{pde, bc}, T[x, y], {x, y}, Assumptions -> {a > 0, b > 0}]


Where DSolve does not solve it, do the following

bc = {T[0, y] == 0, T[a, y] == 0, T[x, 0] == 0, k*Derivative[0, 1][T][x, b] == h*(T[x, b] - T3)}
sol1 = DSolve[{pde, bc}, T[x, y], {x, y},   Assumptions -> {a > 0, b > 0}];

bc = {T[0, y] == 0, T[a, y] == 0, T[x, 0] == T2,k*Derivative[0, 1][T][x, b]  == 0}
sol2 = DSolve[{pde, bc}, T[x, y], {x, y},   Assumptions -> {a > 0, b > 0}]

bc = {T[0, y] == 0, T[a, y] == T2, T[x, 0] == 0, k*Derivative[0, 1][T][x, b]  == 0}
sol3 = DSolve[{pde, bc}, T[x, y], {x, y},   Assumptions -> {a > 0, b > 0}]

bc = {T[0, y] == T1, T[a, y] == 0, T[x, 0] == 0, k*Derivative[0, 1][T][x, b] == 0}
sol4 = DSolve[{pde, bc}, T[x, y], {x, y},   Assumptions -> {a > 0, b > 0}]


And the solution is

sol = First[(T[x, y] /. sol1)
+ (T[x, y] /. sol2)
+ (T[x, y] /. sol3)
+ (T[x, y] /. sol4)]


• Seeing your solution here prompted me to ask a similar question I have been facing recently. I have asked it here. Can you have a look at it. Hope referring in this way is not violating the SE code of conduct. Jul 22, 2020 at 9:52
• FYI, none of your four NDSolve commands (with the correct Robin BCs) return solutions in Mathematica v12.0. I suppose that implies that Wolfram developers have been improving these algorithm recently. Perhaps there's hope for v12.2. Jul 22, 2020 at 15:56
• @MichaelSeifert I suppose that implies that Wolfram developers have been improving these algorithm recently yes, DSolve have been improving. some charts are here is-mathematica-adequately-improving Jul 22, 2020 at 21:43

The primary goal is to solve the pde symbolically.

Then forget separation of variables, try finiteFourierSinTransform. First interpret the PDE and corresponding b.c.s to Mathematica code.

With[{T = T[x, y]}, eq = D[T, x, x] + D[T, y, y] == 0;
bcx = {T == T1 /. x -> 0, T == T2 /. x -> a};
bcy = {T == T2 /. y -> 0, k D[T, y] == h (T - T3) /. y -> b}];


Then introduce the transform $$T(x,y)=u(x,y)+\frac{T_2-T_1}{a}x+T_1$$ to make the b.c.s in $$x$$ direction homogeneous. This isn't actually necessary, but will make the resulting series solution converge faster.

transform = T -> Function[{x, y}, u[x, y] + (T2 - T1)/a x + T1];

{neweq, newbcx, newbcy} = {eq, bcx, bcy} /. transform // Simplify


Eliminate derivative in $$x$$ direction using finiteFourierSinTransform:

tset = finiteFourierSinTransform[{neweq, newbcy}, {x, 0, a}, n] /. Rule @@@ newbcx


Solve the resulting ODE system:

tsol = u[x, y] /.
First@DSolve[tset /. HoldPattern@finiteFourierSinTransform[a_, __] :> a, u[x, y], y] //
Simplify


Transform back:

sol = T[x, y] /. transform /.
u[x, y] -> inverseFiniteFourierSinTransform[tsol, n, {x, 0, a}]


Notice I've used C to denote $$\infty$$ in sol.

To check its validity, compare it to the numeric solution:

xR = 13; yR = 17;

para = {k -> 2, h -> 3, T1 -> 5, T2 -> 7, T3 -> 11, a -> xR, b -> yR};

test[x_, y_] = sol /. para /. C -> 10 // ReleaseHold;

nsol = NDSolveValue[{eq[[1]] == NeumannValue[-h/k (T[x, y] - T3), y == b], bcx,
bcy[[1]]} /. para, T, {x, 0, xR}, {y, 0, yR}]

Manipulate[Plot[{nsol[x, y], test[x, y]}, {x, 0, xR},
PlotStyle -> {Automatic, {Red, Dashed}}, PlotLegends -> {"FEM", "Series(10 terms)"},
PlotRange -> {0, 15}], {y, 0, yR}]