# Solving a 2-D Laplace equation with boundary conditions

I am having a problem solving a 2-D Laplace equation with 4 boundary conditions. I want to solve this differential equation(pde) of phi2[x,y] in Fourier series.

Note that phi1[x] is previously solved 2nd-order polynomial function of x.

phi1[x_] := (eox esi nd q tch tox + esi^2 nd q tox^2 - 2 eox^2 esi vfb + 2 eox^2 esi vg + eox^2 nd q tch x - eox^2 nd q x^2)/(2 eox^2 esi)


To solve pde, I write in my code:

eqn2 = D[phi2[x, y], {x, 2}] + D[phi2[x, y], {y, 2}] == 0
bc3 = phi2[x, 0] == vs - phi1[x]
bc4 = phi2[x, l] == vs + vd - phi1[x]
bc5 = phi2[tch + tox esi/eox, y] == 0
bc6 = phi2[-tox esi/eox, y] == 0

Ω = Rectangle[{-tox esi/eox, 0}, {tch + tox esi/eox, l}]

DSolve[{eqn2, bc3, bc4, bc5, bc6}, phi2[x, y], {x, y} ∈ Ω]


I expect the solution would be in forms of Fourier series. However, the output only repeats my input, not solving differential equation.

MMa is still in its infancy as far as solving pde's analytically. Solving this problem by separation of variables is fairly straightforward, although somewhat messy.

Clear["Global*"]


I will simplify your psi1

phi1[x_] = a + b x + c x^2;


The pde

eqn2 = D[phi2[x, y], {x, 2}] + D[phi2[x, y], {y, 2}] == 0


The boundary conditions

bc3 = phi2[x, 0] == vs - phi1[x];
bc4 = phi2[x, l] == vs + vd - phi1[x];
bc5 = phi2[tch + d, y] == 0;
bc6 = phi2[-d, y] == 0;


Assume

phi2[x_, y_] = X[x] Y[y]

eqn2/phi2[x, y] // Expand
(* X''[x]/X[x] + Y''[y]/Y[y] == 0 *)


If we move the Y part to the rhs we see that the X part is the negative of the Y part. Each side must be equal to a constant. Let the X part be sinusoidal. So we have.

xeq = X''[x]/X[x] == -alpha^2

yeq = Y''[y]/Y[y] == alpha^2


DSolve the xeq for X[x]

DSolve[xeq, X[x], x] // Flatten;
X1 = X[x] /. % /. {C -> c1, C -> c2};


And the yeq

DSolve[yeq, Y[y], y] // Flatten;
Y1 = Y[y] /. % /. {C -> c3, C -> c4};


Convert the above to trig

ExpToTrig[Y1] // Collect[#, {Sinh[alpha y], Cosh[alpha y]}] &
(* (c3 - c4) Sinh[alpha y] + (c3 + c4) Cosh[alpha y] *)

Y1 = % /. {c3 - c4 -> c3, c3 + c4 -> c4}
(* c3 Sinh[alpha y] + c4 Cosh[alpha y] *)


We have the pieces. Put them together as

phi2[x_, y_] = X1 Y1;


To start, look at bc5

bc5
(* (c3 Sinh[alpha y] + c4 Cosh[alpha y]) (c1 Cos[alpha (d + tch)] +
c2 Sin[alpha (d + tch)]) == 0 *)

Solve[bc5, c2]
(* {{c2 -> -c1 Cot[alpha (d + tch)]}} *)

c2 = c2 /. %[];


Now bc6

bc6
(* (c3 Sinh[alpha y] +
c4 Cosh[alpha y]) (c1 Sin[alpha d] Cot[alpha (d + tch)] +
c1 Cos[alpha d]) == 0 *)

bc6[[1, 1]] // Collect[#, c1] &
(* c1 (Sin[alpha d] Cot[alpha (d + tch)] + Cos[alpha d])*)

%/c1;

alphaeq = % == 0
(* Sin[alpha d] Cot[alpha (d + tch)] + Cos[alpha d] == 0 *)


The above determines values of alpha. It is a transcendental equation solved numerically, once constants of the problem are assigned. There are an infinite number of then and the solution becomes an infinite sum over the values of alphas. Ignore the sum for now.

phi2[x, y] // Simplify
(* c1 Csc[alpha (d + tch)] Sin[
alpha (d + tch - x)] (c3 Sinh[alpha y] + c4 Cosh[alpha y]) *)


Combine constants. c1 is redundant

c1 = 1

phi2[x_, y_] = phi2[x, y] // Simplify


Now

bc3
(* c4 Csc[alpha (d + tch)] Sin[alpha (d + tch - x)] == -a - b x - c x^2 +
vs *)


Use orthogonalities of the alphas to solve for c4.

c4*Csc[alpha*(d + tch)]*Integrate[Sin[alpha*(d + tch - x)]^2, {x, -d, tch + d}] ==
Integrate[(-a - b*x - c*x^2 + vs)*Sin[alpha*(d + tch - x)]^2, {x, -d, tch + d}]//Simplify;

c4rule = Solve[%, c4] // Simplify;


Keep thing simple. Don't assign it yet.

bc4

(* Csc[alpha (d + tch)] Sin[
alpha (d + tch - x)] (c3 Sinh[alpha l] + c4 Cosh[alpha l]) == -a -
b x - c x^2 + vd + vs *)


Now orthogonality to solve for c3.

Csc[alpha*(d + tch)]*(c3*Sinh[alpha*l] + c4*Cosh[alpha*l])*
Integrate[Sin[alpha*(d + tch - x)]^2, {x, -d, tch + d}] ==
Integrate[(-a - b*x - c*x^2 + vd + vs)*Sin[alpha*(d + tch - x)]^2, {x, -d, tch + d}]//Simplify;

c3rule = Solve[%, c3] // Simplify;


Now assign c3, c4

c4 = c4 /. c4rule[] // Simplify;

c3 = c3 /. c3rule[] // Simplify;


We have solved for all the constants satisfying the boundary conditions which gives us phi2[x,y] in terms of alpha. The total solution is an infinite sum over all the alphas, determined by numerically solving the equation for alpha` once numerical values for all your constants are specified. I have not printed the output for the later values of these calculations because they are lengthy.